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CHAPTER 1 

INTRODUCTION 

Nickel-cadmium storage catteries have been in use for 
more than fifty years: since World War II, r.ew design princi¬ 
ples have been introduced which have greatly extended their 
usefulness. Many theories and investigations have been put 
forward concerning the reaction mechanisms of both cadmium and 
nickel electrodes. Because the activity of the nickel electrode 
has been more fully characterized than that of the cadmium 
electrode, we will focus our attention on the behavior of the 
cadmium electrode. 

There are a number of methods used for the fabrication 

of cadmium electrodes. These methods are usually classified 

into three categories: (a) pressed powder or pasted type v ^, 

(4 5 6 7 S) 

(b) impregnated nickel sinter type v , - 3t ’ ' ' arid (c) electro- 

(9 10 11 12) 

chemical impregnation 71 ’ ’ The last method, however, is 

simple in operation and can, in theory, achieve high loading of 
active material in a single impregnation cycle at the lowest cost. 

In the electrochemical impregnation process, electro- 
chemically active material is introduced into the porous nickel 
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plaqe by electrochemical precipitation. The plaque, which serves 
as a cathode, is submerged into an electrolysis cell containing 
an aqueous solution of cadmium nitrate at a temperature near 
the boiling point of the solution. The initial pH value of the 
impregnation solution is adjusted to around 3 by titration with 
nitric acid. (As the pH value of the bath increases, i.e., 
becomes more basic, cadmium hydroxide begins to precipitate in 
regions other than within the pores of the plaque, and active 
material is lost in the bath, which is not what we want.) 

The electrochemical reaction at the cathode during impre¬ 
gnation could either be the liberation of hydrogen, 

2H 2 0 + 2e" -► H 2 + 20H" , (1) 

or the reduction of nitrate ion, 


( 2 ) 


3oth reactions produce hydroxy ions which then either diffuse 
away from the electrode porous surface or react with the cadmium 
ions present in the pores to form cadmium hydroxide. The cadmium 


H 2 0 + NO^ + 2e“-» N0 2 " +20H 

• • • • • 

2H £ 0 + NO^" + 2e” - HNC> 2 +30H 

t 9 • • • 

6H 2 0 + ;i0 3 " + 8 e~-- NH +90H 
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hydroxide then deposits within the pores of the porous plaque. 

The precipitation reaction is: 

Cd +2 + 20H“ -* Cd( OH) 2 . (3) 

There is the possibility of forming cadmium according to the 
reduction reaction: 

Cd(0H) 2 + 2e~ -*- Cd + 20H". (4) 

This reaction is presumably prevented from taking place by 
hydroxy ion generated on the electrode surface by reaction (1) 
or (2), so that there is no chance for cadmium ions to diffuse 
through the solution onto the surface and be reduced to cadmium 
metal. 


Following impregnation, the electrode is cathodically 
polarized in a 25 $ potassium hydroxide solution with the impre¬ 
gnated plate connected as the cathode and subjected to low direct 

2 

anodic current (0.25 to 0.75 amp/in ). The solution is maintained 
at about 100°C for a period of from 15 to 45 minutes to convert 
cadmium hydroxide to cadmium in the pores of the plaque by 
reaction (4) or to remove any nitrate ion or other undesirable 
reduction products. The active material of the electrode is pre¬ 
sumed to be both cadmium hydroxide and cadmium. 
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In this work, the electrochemical impregnation process 
for the fabrication of the cadmium electrode was examined both 
experimentally and theoretically. 

Experimentally, a set of current-density-vs.-time transients 
at constant applied potential were obtained for a nickel micro¬ 
electrode. The current density was due to the occurrence of the 
electrode reaction (2) and precipitation reaction (3) presented 
earlier. For the purpose of identifying the electrochemical 
reactions, another set of current-density-vs.-time transients were 
obtained using an electrolyte which contained no cadmium ion. 

In this case, the precipitation reaction (3) no longer took 
place. The reactions at the electrode were the only reactions 
taking place. 

The transient current responsed in an electrolyte which 
contained cadmium ions differed significantly from that in an 
electrolyte which did not contain the cadmium ions. This diff¬ 
erence is presumed to be due to the precipitation reaction. 

The current obtained in the presence of cadmium ions in the 
electrolyte at any time was higher than that obtained at the 
same applied potential in the absence of cadmium ions. This 
result is explainable by the fact that the cadmium ions in the 
solution consumed the hydroxy ions which were produced by the 
heterogeneous reaction (2) and, as a result, promoted the hetero¬ 
geneous electrode reactions. 
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The experimental current-density- vs.-time curves obtained 
in the potential-step experiment were used to identify both 
the heterogeneous electrochemical and the homogeneous preci¬ 
pitation reaction rate expressions. The heterogeneous rate was 
identified first from the experimental results in which there 
was no cadmium ion in the electrolyte. The precipitation rate 
was then identified from the experimental results in which the 
cadmium ions were present in the electrolyte. 

A set of differential equations which described the trans¬ 
port process in the cadmium-ion-free solution was solved in terms 
of two unknown constants. The unknown constants, related to the 
heterogeneous electrode reactions, were identified by matching 
the theoretical results with the results of the experiment. The 
precipitation reaction rate constant was then obtained by a 
similar procedure. 
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CHAPTER 2 


LITERATURE REVIEW 

Impregnation of nickel sinters with cadmium hydroxide, 
followed by electrochemical reduction to cadmium, is the most 
economical method by which to fabricate cadmium electrodes. 

Almost all of the research works concerning this subject may 
be divided into three typesi (a) manufacturing techniques de¬ 
velopment, (b) impregnation process study arid (c) charge-discharge 
study. 


During recent years, manufacturing techniques, especially 
those involving electrochemical impregnation, have employed 
widely varying operating conditions^. There are 
several parameters that affect the electrochemical impregnation 
process such as: 

1) Temperature: Beauchamp emphasized the fact that high 

impregnation bath temperature keeps the size of cadmium 
hydroxide crystals small, and thereby obtaining a more active 
electrode. 

2) pH value: In Beauchamp’s process , the pH value was 

closely controlled by buffering the cadmium nitrate solution 

/ A ' \ 

with sodium nitrite. Pickett °' in a similar process, 
maintained the pH value of the cadmium nitrate bath between 
3-5 by titrating the solution with proper amount; of nitric 
acid. 



7 


3) Current density: The current density, based on apparent 

(14) 

electrode surface, is another important variable. Beauchamp 
used D. C. current of 4-8 amps/in^ while Bulan ^-5) usec j 

2 

alternating current of 1.2-1.6 amps/in and claims that the 
electrode impregnated by alterning current technique retained 
more of its capacity after cycling as compared to electrodes 
fabricated by other methods. Pickettp resen ts another 
technique similar to the above and claims that his technique 
yields higher loading cadmium electrodes and achieves about 
85% utilization of the active material. 

The electrochemistry of the Cd/Cd(0H) 2 electrode in con¬ 
centrated alkaline solution, typically KOH as in the Ni-Ca 
battery, has been studied extensively both theoretically and 
experimentally ^^. Bennion, et al.^^ in their 

theoretical development, present two models, the solution- 
diffusion model and film model, in considering the failure 
mechanism in the performance of secondary batteries using metal/ 
metal salt couples. They conclude that one mode of failure of 
the metal/metal oxide porous electrode is the blocking of the 

pore or the complete covering of surface by product crystallities. 

( 18 ) 

Falk' obtained X-ray diffraction patterns from electrodes 
submerged in electrolyte during charge and discharge by means 
of a special test cell. He claims that, during charge, Cd(CH) 2 

is successfully transformed into Ca metal in the cadmium electrode, 
but even after a strong overcharge, this transformation is not 
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complete, as indicated by the X-ray pattern at this stage. 

Furthermore, Cd(OH), is the end product during discharge and 

there is no evidence of the presence of a CdO film which is 

believed to be responsible for the passivation of the electrode. 

( 19 ) 

Gross and Glocking summarized the results of recent research 
by characterizing the negative cadmium electrode in a nickel- 
cadmium cell and evaluating some of the published results about 
which researchers still disagree among themselves, such as the 
mechanism for passivation, which is an important phenomena that 
affects the capacity of the nickel-cadmium cell. 

In comparison to the extensive studies on the electro¬ 
chemistry of the Cd/Cd(OH) 2 electrode in KOH solution, there 
is little information available concerning the electrochemical 
reactions taking place during the electrochemical impregnation 
process which uses cadmium nitrate solution as electrolyte. 

The only research works known to the author are the electro- 

(21 22 ) 

cnemical impregnation studies using cycling voltammetry ’ 

2 2 
They used a 0.013 cm nickel microelectrode and a 0.064 cm 

cadmium microelectrode. Some important conclusions of their 

studies are summarized below: 

(1) There was a large reduction wave at -O. 63 V vs SGE(Saturated 
Calomel Electrode) which appeared only during the initial 
negative scan. This is believed to be the reduction of a 

solution species to liberate 0H~ ion in the presence of 

+2 

Cd to precipitate cadmium hydroxide on the electrode 


surface. 
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(2) A second reduction process was observed at -0.62V vs SGE 
during the reverse positive scan. This is due to the reductive 
formation of metallic cadmium from cadmium oxide or hydroxide 
that was deposited on the electrode surface during the pre¬ 
vious forward scan. 

(3) Some processes (or one of the processes) during the first 
cycle effectively passivated the electrode to further electro¬ 
activity. The formation of the gray material at -O. 63 V vs 
SCE, discussed in (1), was the most probable reason for the 
passivation of the electrode. 
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CHAPTER 3 

THEORETICAL APPROACH 

3.1 General 

In voltammetry at constant voltage, the current through 
the electrolytic cell is measured as a function of time while 
the potential of the polarized electrode is held at a constant 
value. The current is affected by both the transfer process and 
the electrochemical and/or chemical reactions. 

A plane electrode is immersed in a solution which contains 
a substance 0. At the imposed potential, 0 is reduced to another 
substance R at the electrode surface, according to the following 
electrochemical reaction: 

0 + ne" -► R (5) 

The current-density-vs.-time curve obtained at a constant 
voltage depends on both the kinetics of the electrochemical 
reaction and the rate of mass transfer. 

Under our experimental conditions, diffusion is the only 
mode of mass transfer. The other two modes of mass transfer, 
namely, migration and convection, are insignificant. Migration 

I 
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is eliminated by the presence of a large excess of supporting 
electrolyte which reduces the potential field. Convection is 
avoided due to the conditions that (1) the solution is not 
stirred and (2) the duration of electrolysis is short, so that 
the density change has not yet become a significant factor and 
natural convection has not taken place. 


The dimension of the cell, as compared to the microelectrode 
used, is so large that the solution may be regarded as extending 
to infinity; i.e., the diffusion process is semi-infinite. 
According to Fick's first law, the flux of the substance 0 at 
a distance x from the electrode and in the direction perpen¬ 
dicular to the electrode is: 

^C 0 (x,t) 

N 0 (x.t) = - D 0 -“ (6) 

OX • 

Conservation of the species 0 leads to the following differential 
equation: 

?>C n (x t t) = D ^ 2c 0 (x,t) 

t 0 7) x 2 (?) 


where Cq is the concentration of substance 0, which is a function 
of x and t. Dq is the diffusion coefficient of substance 0, 
and is assumed to be a constant. 
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When the substance 0 is reduced at a very high rate, by 
imposing a sufficiently large potential step, such that the 
concentration of this substance at the electrode surface is 
reduced to zero soon after the start of electrolysis, the 
boundary condition at the electrode/electrolyte interface for 
equation (?) is CQ(0,t)=0 for t>0. Initially before the electro¬ 
lysis, the solution is homogeneous, and the concentration Cq 
at t=0 is uniform with the bulk value. The other boundary 
condition is that C Q approaches its bulk value when x is 
sufficiently large. 

The electrolysis current is equal to the product of the 
charge involved in the reduction of one mole of substance 0 by 
the flux of this substance at the electrode surface. Thus the 
current density (current per unit electrode surface area) is: 

1= nFN 0 (0,t), (8) 

where n is the number of electrons involved in the reduction 
reaction .of substance 0, F is the Faraday constant, and NQ(0,t) 
is the flux of substance 0 for x=0. 

The solution of the above problem is^-^ : 

1= nFD * cj -x - ~ ”- (9) 

U U t 2 

f 

where Cq is the bulk value of substance 0. 
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The current density represented by equation (9) is called 
the limiting current density of electrochemical reaction (5), 

This is the maximum current density which can be achieved by 
this reaction. It was found that if nitrate ion is the only 
reducible substance, the current density predicted by equation 
(9) is considerably lower than the experimental value. This 
contradition is resolved if one realizes that under the experi¬ 
mental conditions water is also reduced. 

At any other less-negative potential the current will 
depend on the electrode kinetics. Now the condition thax the 
reducible substance concentration becomes zero after the appli¬ 
cation of the potential step has to be replaced by a kinetics 
equation, which is a function of potential and the concentrations 
of surface reactant and product. No general solution is available 

except for the case when the reaction rate is linearly dependent 

, . ( 21 ) 
on tne concentrations of the surface reactant and product J . 

3.2 Description of the Model 

The real porous electrode is formed by sintering nickel 
powder on a nickel grid to form a plaque. The plaque serves as 
the cathode, which is positioned in a cadmium nitrate solution 
with a pH value between 3 to 5* The passage of the current 
produces hydroxy ions which cause the cadmium to be deposited 
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into the porous pores. The transport phenomena associated 
with electrochemical reactions in a porous medium is very 
complicated. The approach here is to study the intrinsic 
reactions which take place inside the pores in a simple 
flat electrode. The plate electrode is submerged in a 
semi-infinite pool of electrolyte so that the mass transfer 
problem can be treated as a one-dimensional problem in 
the x-direction (perpendicular to the electrode surface). 

In the solution side, there is a thin double-charge layer 
near the electrode. The thickness of the double layer is 

O 

on the order of 10 to 100 A. It is too thin to be probed 
adequately, and the theory of the diffuse layer is a 
microscopic, rather than a macroscopic, phenomena. The 
mass transfer region to be considered here is thus located 
outside of this double layer. 

The electrodeposition process operates on the principle 
that the hydroxy ion, which is generated by the reduction 
process, coprecipitates with the cadmium ion in the solution 
to form cadmium hydroxide crystalloid. The approach taken 
here to unravel the complicated sequence of reactions is 
to form a model which describes the transport processes 
and reactions. The model will first be used to analyze 
the experimental results obtained under the condition of 
no pricipitation reaction, by using an electrolytic solution 
which is free of cadmium ions. The heterogeneous reduction 
reaction will be determined. The transport and reacxion 
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problems will then be examined in the case when the solution 
contains cadmium ions. The rate of the precipitation 
reaction will then be determined by matching the model's 
prediction and the experimental results. 


3.3 Formulation of Transport Process Equations 


Three ionic species to be considered in the electro- 

_ _ +2 . 

chemical system are NO^, GH and Cd ions. In the absence 
of migration and convection, the flux of each species is 
expressed as^^^s 


N i - " D i 


, 1 = 1 ,< 2,3 , 


where 


is the flux of species i (moles/cm -sec), 
is the diffusion coefficient of species i (cm /sec), 
Ch is the concentration of species i (moles/cur ), 

which is a function of both x and t, and i=l,2,3 

- - +2 

corresponding to NO^, GK and Cd ions. 

(24) 

The equation of continuity of each species is 


? C i = 

"St 


d N i 


i=l,2,3 , 


where rn is the rate of production of species i 
(raoles/cm -sec) by precipitation reaction, which occurs 
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only in the presence of cadmium ions. 


Substitution of equation (10) into (11) give: 


?>t 


2 

3 c. 


= Di' 


2 ):. 


2 


i=l,2,3 


( 12 ) 


It is assumed that the precipitation reaction ( 3 ) 
can be treated as a homogeneous reaction and is linearly 

proportional to the degree of supersaturation. Therefore, 

the consumption rate of cadmium ion due to this reaction 

is ; 


R 


3 



(13) 


where k is the homogeneous reaction rate constant (sec - ^), 
Ksp is the solubility constant of Cd(0H ) 2 

and has the value of 2 xl 0 - ^ (moles^/cm^) ^5 ). 

3y the stoichiometric relation of equation (3)> the con¬ 
sumption rate of hydroxy ion due to the same reaction is : 


«2 = 2 R3 


-2k 




(c 2 )‘ 


(14) 




Substitution of equations (13) ana (14) into equation (12' 
leads to the following two differential equations s 


7)t 


->c. 


- d ° 2 = 2 k 


Ksp 




(c 2 )‘ 


( If) 



2 { 

T'3 

D-.-k C_ 

3}x 2 I 3 


t p \ 2 
v 2 ; 


For NO" ion, 


R 1 = o. 


Equation (12) for NO^ ion may be reduced to 


2> C 1 2T c i 

-— = n -— 

■at i \ 2 

O X 


Only two of the three differential equations need 

to be solved for the concentration profiles, due to the 

(24) 

electroneutrality condition s 


2C 3 -C 1 -C 2 =0 


Substitution of equation (18) into equation (15) 


gives s 


~2> C 2 
it 






Equations (1?) and (19) are solved simultaneously 

to obtain the concentration profiles of ar.d C~, i.e., 

- . +2 

i.'O^ and OH ions. The profile can then be obtained by 
equation (18) from the known profile^ of and CU. 


A set of differential equations, which are simpler 
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than equations (17) and ( 19 ). ars obtained when the elec 
trolyte contains no cadmium ion. In the absence of cad¬ 
mium ion, the precipitation reaction ( 3 ) can not occur; 
thus, FL terms are all equal to zero. Equation (12) is 
then reduced toi 


7> C i 

at 



i=li2 , 


( 20 ) 


in which i=l,2 corresponds 


to N 0 ^ 


and OH 


ions. 


Equation (20) is written for and C^s 


D C 1 

7>t 



= 

l>t 



( 21 ) 


( 22 ) 


The other condition, the concentrations of various 
ions at the electrode/electrolyte interface after the 
start of the current flow, has to be described by the 
electrode kinetics of the reduction reaction of N 0 “ ions 


3.^ Electrode Kinetics 


The possible reduction reactions are numerous, as 
shown in equation (2). Since in e*..ch reaction, OH" ion 
is generated, we will represent the production of OH" 
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ions as the following electrochemical reaction: 

2H 2 0 + N0^ T 2e“- * HN0 2 + 30H". (23) 

Because there are two electrons involved in the reaction, 
it is likely not an elementary step. Let's assume the 
reaction comprises two elementary steps, each step involving 
the transfer of one electron, and a dissociation reaction: 



(slow), (24a) 

(fast), (24b) 

(fast), (25) 


where and K c2 are cathodic reaction constants, and K ^ 
and K a2 are anodic reaction constants, respectively. is 
the dissociation reaction constant. (OH ) 2 ion is an unstable 
intermediate ion which immediately enters the next reaction (24b/ 
The rates of these elementary steps should always be propor¬ 
tional to each other. For the elementary steps listed above, 
reaction (24b) should occur once every time reaction (22a) occurs 


Let I. and I, denote the cathodic current densities 
of reactions (24a) and (24b), and r, and r 2 denote the 
reaction rates of reactions (24a) and (24'c), respectively. 
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Then, tne relationship between r. , I- and surface ion 


concentrations of N cZ, (CK) 2 , which app 

(2 ^) 

step (24a), may be written as 


ear in e-emer.tar; 


- _Ll 

'l" nF 


(n=l; 


3 K al ( C NC;]( C (0H)-) f J 

W W 


(:-/3/:? 


.. r„o I f r o 1 o „^ r f 

cl °ho: C H 9 0 e '‘ ? " " 

u JJ u ^ L *" 


r a? „ 


Similarly, for reaction (24b), one has: 


(n=l) 


^K-r«pp^v] 

K o 2 [ C (0K)-] ex? [-^ Lv j' 


where (3, and ^ are symmetry factors representing the fractions 
of applied potential V, which promotes the cathodic reaction. 
C.°.~- , and Z%,- are the surface ion concentrations of 

4'0y (OH)^ and OH” on the electrode surface. 


Since reactions (24a) and (2^b) occur at the same rate, 
we have: 


1= 1^ f I, s 2I 2 = 21. , 


2 
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where I is overall current density. It is further assumed 
that reaction (24b) is fast, and is essentially in equilibrium. 
From equation (2?), 


-r 

x c2 ex p ['-fr 7 ] 


and = 
r 


=r = overall reaction rate. 


Substituting equation (28) into equation (25) and 
combining equation ( 29 ) gives: 


_ K al K a2 ( r o 


[ C °°-][ C0 OH-] 2 exp £ 


(2-A)? 


x d( c kJfS 2 o]-p[-^ v ] . 


From reaction (25) one may write 


K d = 


( G HN0 2 ][ g 0H~] 

[ c no“][ c h 2 oJ 


Substitution cf the above expression into equation ijC) leads 
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ions and electrolyte containing no cadmium ions, two 
and four boundary conditions are needed to solve the 
equations (17) and (19) or (21) and (22). 


initial 
differen 


At t=0, the concentration of each ion is equal to its 
bulk value at each point along the x-direction. The initial 
conditions are as follows: 

1# at t=0, C^= for all x, 

2. at t=0, C 2 = C 2 for all x, 

where the superscript "c" denotes the bulk condition. 


The first boundary condition is the consequence of 
the fact that the flux at the interface has to be equal to 
the rate of heterogeneous electrode reaction, i.e.: 

xm0 * K 2 C c °) 3 - X lC°?l <32) 

where superscript "o" denotes the surface condition. To 
satisfy the stoichiometric relation of equation ( 23 ), one 
has, at x=0: 


* ^ 

, > : 1 

- - 


f \ 

> C 2 

dx 

1 

o 

II 

X 

^2 -v 

3x 

^ J 


k J 


This is the second boundary condition. The ct.-.er two 
boundary conditions are: 
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3 • at x=oo , C^= C°, 

•u 

4. at x=co , C 2 = C 2 . 

That is, the concentrations far away from the electrode/ 
electrolyte interface do not change. 


The flux in equation (32) is proportional to the curren 
density as a function of time: 
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CHAPTER lr 

EXPERIMENTAL WORKS 

A potential-step was generated from a function 
generator and was applied to the working electrode/elec- 
trolyte interface in the electrochemical cell. The ca¬ 
thodic current was thus sampled and stored in memory at 
a constant time interval by a microcomputer. The resulting 
current/time transient was then plotted by a x-y recor¬ 
der or printed out by a line printer by recalling these 
current data from the memory of the microcomputer. 

4.1 The Electrochemical Cell 


The current transient under a potential step con¬ 
dition was obtained in a three-electrode cell. The cell 
and electrodes configuration is sho-wn in Figure 1. The 
working electrode was a nickel microelectrode. It was made 

of a nickel wire which was pressure-fitted into a teflon 

. 2 
cylinder. The flat, exposed surface was 0.013 cm 4 ". A 

saturated calomel electrode (SCE) and a cadmium bar were 
used as the reference and che counter electrodes, res¬ 
pectively . 


The cell was filled with electrolyte of either 
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cadmium nitrate solution or potassium nitrate solution. 

The pK value of the solution was controlled. Potassium 
chloride was used as the supporting electrolyte. The electro¬ 
lytic solution was not stirred during the experiment, so that 
the conditions of semi-infinite linear diffusion were maintained. 

Two experimental parameters were varied: 

1) the ion concentration of the solutions and 

2) the magnitude of the potential-step applied to the 
interface. 

4.2 The Instruments and Electrical Setup 

The-block diagram for the experiment setup is shewn in 
Figure 2. The instruments shown in Figure 2 contain the 
following equipment: 

a) A Princeton Applied Research (PAR) Model 1'~’3 Potentiostat/ 
Galvanostat : In the experiment, the operating mode was set 
at CONTROL E., which means that the instrument functioned 
as a potentiostat. The current was measured while keeping 
the potential of the working electrode constant. This was 
accomplished by setting the counter electrode to the 
required level so that the 'working electrode potential was 
at a programmed potential with respect to the reference 
electrode. The instrument is provided with a Me del l"o 
Current-to-hoitags Converter which is the "basic" piug-in 
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module. The M 176 provides a do voltage output ’./hioh is 
proportional to the current. The cell current can be eitner 
displayed-on the panel, by a x-y recorder, or monitored by 
the microcomputer. Connection to the cell was made through 
the external cable» The counter electrode was connected to 
the red clip, the working electrode to the green clip and 
the reference electrode to the Electrometer Probe which has 
a very high impedance, thereby insuring us that the current 
was flowing to the reference electrode, 

b) A PAR Model 175 Function Generator ! This is a programmable 
waveform generator which has two operating modes, SWEEP and 
PULSE. The former generates'a sequence of triangular waves, 
while the latter generates a sequence of step function. Eor 
the potential step experiment, the operating mode was set-a 
PULSE mode and the pulse width selector at STEP position. 
The magnitude of the step-function was set by the setting o 
the 3 potential (upper limit) in the POTENTIAL section on 
the front panel of this instrument. 

c) An Intel CPU S080A Bases Microcomputer (cU-K ?,A The micro 
computer was the central part of the experimental setup. 

The functions of the microcomputer in the experiment were: 

(1) to trigger the Function - 'Generator which activated the 
potential step presented earlier, 

(2) ...to sample the current outputs from the M 17c Converter 

- -and store-them in-the-memory. — 

(3) to convert these digital data back into analog signals 
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and then plot them on the x-y recorder. 

(4) to print the digital data on the line printer. 

(5) to store these digital data in the tapes. 

The peripheral equipment of the microcomputer include a 
CRT (TV screen), the keyboard, a cassette recorder and a 
line printer. 

d) A Houston Model RE 0074 x-y Recorder : The Recorder receives 
the analog outputs from the DAC (Digital-Analog Converter) 
and plots them on the graph paper. 

4.3 The Computer Program 

A computer program ’written in BASIC language -was used no 
carry out the experiment. The program consists of a main program 
and a machine-language subroutine called "MUG". This subroutine 
includes an execution statement which generates a trigger signal 
to activate the Function Generator. A potential step was then 
immediately applied to the working electrode, and thereafter 
the computer started to sample and store the current output at 
a fixed time interval. The analog signal was converted into a 
digital datum by a 12 bits ADC (Analog-Digital Converter) and 
then stored in the memory. When the number of samples reached 
a preset value, the sampling routine was terminated. The digital 
sigr.als stored in the memory were binary numbers. These data 
were converted into decimal numbers and then converted into 
analog signals,, and finally plotted on the x-y recorder and 
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printed on the line printer. This was executed by the main 
program. The main program and the "MUG" subroutine are 
presented fully in Appendix A1 and A2. 

4.4 Experimental Conditions and Procedures 


In the experiment, the initial pH value of the electro¬ 
lyte was maintained at 3*G, which was the condition used in the 
actual electrochemical impregnation process. This was 
accomplished by titrating the electrolytic solution with a 
concentrated KC1 solution. The concentrations of nitrate 
ion used were 0.005M, 0.002574 and 0.00125M prepared from either 
the KNC^ or Cd(NO^) solution. The cadmium ion concentrations 
were 0.0025M, 0.00125M and 0.00Q625M, All the electrolytic 
solutions contained 0.5M KCi as the supporting electrolyte. 

It was .found later that the electrolysis of water made 
significant contribution to the total current, because the 
current obtained was significantly higher than the limiting 
current. In order to obtain the current due only to the reduction 
of nitrate ions, several .runs with solutions which contained no . 
nitrate ion were conducted. The solutions used for this purpose 
were either water or C.0G25M CdCl^ solution. Beth solutions 
contained 0.5M KCi as the supporting electrolyte. This current 
was subtracted from the current obtained in the oreser.ee of 
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nitrate ions in the solution. The resulting current is the 
faradic current of the reduction of nitrate ions. 

The magnitudes of the potential-step used in this 
study were -0740V, -0.60V and -0.80V versus the equilibrium 
potential which was at approximately -G.35^V vs SCE, 

The experimental procedures are listed blows 

(1) The electrical circuit was set up as shown in Figure 2, 
except that the cell was disconnected from the Potentiostat/ 
Galvanostat by setting the cell selector on the instrument 
to the OFF position. A proper current range was then 
chosen from the Current Range Switch to make sure that 

the current output would not be overloaded. The Function 
Generator was initialized by depressing the INITIAL 
control pushbuttom. The initial potential was set at 
zero volts. 

(2) A known volume of electrolytic solution and an equal 
volume of supporting electrolyte, KC1, were added to the 
cell and titrated with HC1 solution to a pH value of 3*0. 

The solution was then deaerated by bubbling purified 
nitrogen through the stirred solution for about ten 
minutes. The gas continued to pass above the solution 
during one experiment. 

(3) After the working electrode surface was polisr.sd by a 

3/0 alumina paper, it was rinsed thoroughly with distilled 
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and deionized water, and then was positioned in the cell. 

The electrodes (the working electrode, counter electrode 
and reference electrode) were then connected to the 
Potentiostat/Galvanostat as described in section 4.3• 

After the electrodes were correctly connected, the cell 
selector was switched to the EXTERNAL CELL position. The 
rest potential of the solution could be determined by 

turning the knob in the APPLIED POTENTIAL/CURRENT section 
of the front-panel of the Potentiostat/Galvanostat so 
that the current reading displyed on the front panel 
meter was zero. The applied potential was then the rest 
potential. 

(4) The upper limit (B potential) on the Function Generator 
was set. The value of the -3 potential plus the rest 
potential was the total potential applied to the electror 
chemical cell. 

(5) The experiment was initiated by running the computer 
program. The' current/time data were stored in the computer’s 
memory and -then the results were converted to line 

printer output. 

4.5 Experimental Results 

Two sets, a total of twenty-two runs, were made. The 
first set of' runs was conducted in solutions which contained 
no cadmium ion. The second set of runs was made in solutions 
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which contained cadmium ions. For each set of runs, clank runs 
were made in the solutions without the presence of nitrate 
ion for the purpose of determining the background current 
aue'to the water decomposition reaction. The solution for 
blank runs for-the first set was.water, while the solution 
for blank runs for the second set was made by using the 
0 . 0025 M CdC^ solution. 


The current-density-vs.-time data printed by the H 14 
Line Printer are shown in Table 3.1 through Table B.22 in 
appendix 3. The measured current was converted to the current 
density which appears in those tables. This is calculated 
by dividing the current by the electrode surface area 
0 . 013 cm 2 . 


Tables 3.1 through 3.9 are the first data set which used 
potassium nitrate as the electrolytic solution, while Tables 
3.10 through 3,12 are the background currents of water used 
to subtract from the first data set. 


Tables 3.13 through 3.21 are the second data set which 
used cadmium nitrate as the electrolytic solution, while 
Table 3.22' is the background current of water in tr.e presence 
of cadmium ions. 


A series of current-density-vs.-time curves with various 
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electrolytic solution concentrations and applied overpotentials 
are shown in Figures 3 through 6. In these figures, the 
current due to the electrolysis of water has been elirainited 
by using the routine tabulated in Table 1. 

Figures 3 through 5 are the current-density-vs.-time 
curves obtained in the solutions which did not contain cadmium 
ion. Figure 3 shows the current-density-vs.-time curves 
obtained in a solution containing 0 .Q 05 M nitrate ion and at 
potentials -0.40V, -0.60V and -0.80V from the equilibrium 
potential. Figure 4 shows the current-der.sity-vs,-time curves 
obtained in a solution containing 0 . 0025 M nitrate ion, and at 
potentials -0.40V, -0.60V and -0.807 from the equilibrium 
potential. Figure 5 shows the current-density-vs.-time curves 
obtained in a solution containing 0.00125M nitrate ion ax -0.40V, 
-0.60V and -0„80V_from the equilibrium potential. Due to the 
lack of blank runs with O.OOI 25 M and 0.000625M cadmium ion con¬ 
centrations and - 0 . 60 V and -0.80V from the equilibrium potential, 
only one current-density-vs.-time curve was obtained with 0.0G5M 
nitrate ion and 0 . 0025 M cadmium ion concentration at the poten¬ 
tial of -0.40V from the equilibrium potential. This is shown 
in Figure 6. 

4.6 Discussion of Experimental Results 


Some characteristics can be seen from the curves shown 
in Figures 3 through 6. 
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Table Is Evaluation of Currer.t-Density-vs. -Time lata 
Plotted in Figures 3 Through 6. 


— 

Figure is 
plotted 

by subtracting the 
corresponding current 
density at each time 

uoint in 

from that 

of in 


(1) 

. 

Table 3.1 

Table 3.10 

Figure 3 

(2) 

Table 3.2 

Table 3.'11 


(3) 

Table 3.3 

Table 3.12 


(1) 

Table 3.4 

Table 3.10 

Figure 4 

(2) 

Table 3.5 

Table 3.11 


(3) 

Table 3.6 

Table B.12 


i 


( 1 ) 


Table 3.7 


Table 3.10 


























Firstly, the current density decays with tine in a ] 

hyperbolic fashion. This is the result of depletion of the j 

iN'O^ ions near the electrode surface. Theoretically, there j 

is an initial sharp rise of current associated with double-layer 
( 27 ) 

charging , due to the application of a potential-step. This was \ 

not observed in the experiment, however, because the tine for j 

charging is usually very short as compared to the electrolysis 1 

time. 

;i 

i! 

Secondly, the sharp current drops in the beginning indicate j 

that very sharp concentration gradients are established for the i 

: I 

nitrate ions and the hydroxy ions as soon as the surface 

reaction (23) takes place. The concentration gradients are ! 

i 

gradually decreased due to the diffusion layer extending in ^ 

the direction of decreasing concentration gradients. The i 

result is that the decay of the current density is not as fast j 

as in the beginning. 



The curves for the case of potential at -0.80V from the 
equilibrium potential and solutions containing 0 . 00251(1 and 
0.00125M nitrate ion, as shown in the upper-most curves in Figures - 
and 5» respectively, are somewhat different from tne others. 

These two curves do not follow the same trends as their lower 
potential counterparts. It is suggested that the reactant, 
nitrate ion, near the electrode surface, is used up in a very 
short time and another reaction must be taxing place. This 






behavior is the mos* pronounced as shown in the upper-most 
curve in Figure 5* A maximum current is observed at time 
equal to about 0.5 msec. This curve thus strongly supports 
our explanation that some other reaction is taking place at 
that time. On the other hand, the absence of this behavior 
in the curve for the case of 0 . 005 M nitrate ion shown in the 
upper-most curve in Figure 3 suggests that the nixrate ion 
concentration near the electrode surface is not zero. 

Two parameters were varied in the experiment, i.e., 
the concentration of the electrolyte solutions and the value 
of the potential from the equilibrium potential (rest potential). 
The changing of either the bulk concentration of the electrolyte 
solutions or the overpotential that is applied to the working 
electrode will change the value of the current density. For 
a given bulk concentration, the increase in the applied over- 
potential increases the electrochemical reaction rate on the 
electrode surface. The effect of the applied overpotential on 
the current density can also be observed in Figures 3 through 5> 
In each of these figures, a plot of the limiting current - 
density-vs.-time curve is also shown. The. limiting current 


density 

is calcUJ.axed 

by equation (9) 

the nitraxe ion. 

with the 

corresponding 

bulk concentration of 

FecalI 

una 0 xn 0 1 m 1 

current 

is defined as 

maximum currer.x 

one can 0 

ctair* as zte 

applied 

overpoter.xial 

is increased. T 

he surfac 

e c on c en xr axi or. 


of the reactant for the case of limiting current is, in fact, 
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zero at any tine. After examining Figures 3 through 5» find 
that, for the case of -0.30V, some of the values of current 
density exceed the limiting current density. This phenomenon 
becomes more significant when the bulk concentration of pota¬ 
ssium nitrate is decreased. The reason for this phenomenon 
probably is that at the high applied overpotential, many 
reactions can subsequently be expected to take place. The 
maximum current appearing in the curve for the case of -0.50V 
is a result of this behavior. 

When the applied overpotential is fixed, the hetero¬ 
geneous reaction rate constants in both the cathodic and anodic 
directions, K, and K 2 , are fixed. The supply of the nitrate 
ions on the electrode surface comes only in the way of diffusion 
of the nitrate ions from the bulk solution. Increasing the 
bulk concentration will increase the rate of nitrate ion supply 
to the electrode surface, thusly increasing the current density. 
Figures 7 through 9 show the effect of the bulk concentration 
of the electrolytic solution on the current density with fixed 
applied overpotential. From the three figures, we can see that 
the current density is approximately proportional to the bulk 
concentration of potassium nitrate, except for the case of -C : . 

When electrolytic solution contains cadmium ions, a 
different set of results was obtained. The experimental 
current-density-vs.-time curve obtained in a J.C 02pl-l 
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Effect of the Bulk Concentration of the Electrolytic Solution 
on the Current Density with Fixed Applied Overpotential. 
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Cd(N0^)2 solution and at -0.40V applied overpotential is shown 
in Figure 6. In addition to the electrochemical reduction of 
the nitrate ions on the electrode surface, a precipitation 
reaction also takes place in the solution because the cadmium 
ions in the electrolytic solution can coprecipitate with the 
hydroxy ions produced by the electrochemical reaction on the 
electrode surface. The current density is expected to be 
greater than that when the cadmium ion is absent from the 
electrolytic solution. Figure 10 shows the current-density- 
vs.-time curves obtained in a 0.005WI KNO^ solution and in a 
0.0025M Gd(NC ^)2 solution. The applied overpotential is -0.2-07 
for both cases. From the curves, one can see that the current 
density obtained in the 0 . 0025 M Cd(N0^) 2 solution is larger 
than that obtained in a 0.005M KNC^ solution. The reason for 
this is that the consumption of the hydroxy ions due to the 
precipitation reaction establishes a driving force to produce 
more hydroxy ions and thus promote the electrochemical reaction. 
Some interesting facts may be seen from the figure. Firstly, 
the difference between two current densities increases as time 
increases, finally reaching a constant difference; secondly, 
the rate of decrease of the current density with time is slower 
for the case of using Cd(N0^) 2 as the electrolytic solution. 

One may thus obtain a larger steady state current density using 
Cd(N0^ P as the electrolytic solution. 
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CHAPTER 5 

NUMERICAL SOLUTIONS 
OF THE THEORETICAL MODEL 

5•1 Orthogonal Collocation Method 

There are many numerical techniques which can be used 
to solve the differential equations of the theoretical model 
presented earlier in Chapter 3* 

Orthogonal collocation method was selected as the method 
of solution. The orthogonal collocation method requires very 
little implicit mathematics and results in a large savings in 
computer time over the common finite difference scheme such as 
the Crank-Nicholson method. A brief description of the ortho¬ 
gonal coxlocation method is presented in Appendix C. 

5*2 Working Equations 


5*2.1 Dimensionless Form of the Differential Equations 


The orthogonal collocation method requires the region of 
solution to be restricted to the interval C3,13 . This is 
accomplished by normalizing the spatial variable by the 
parameter, £, the diffusion thicxr.ess, ^ =x// . 







dimensionless forms* 

* * ■'-> 
r* _ — , n — - 

1 cl 2 c? 

1 1 

the dimensionless time is defined as: 


Substituting the above dimensionless variables into equations 
(17) and ( 19 ) and rearranging leads to the following dimension¬ 
less eauations: 


.(.Eulh 
l V ^ 2 


* 2„ * 

^ ^ '‘l , , _ * „ *q 3?Ist 

•2T ‘s'? 2 ' r 1 T '“ 2 ; "777r 


where K= -i-£— and 2Ksp= - AS F , 

3 2 

are the dimensionless homogeneous reaction rate constant and 
solubility product, respectively. Similarly, equations (21) 
and (22) become: 


U 1 \ 1> U 1 


* ? * 

z : 2 _ y c 2 


V>J 
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The initial conditions are: 


1. att*0, C. =1, 


for all 0 


s l - 


2. atr« 


=0, C 2 "= C°/C° , for all 0 < £ < 


1. 


The first two boundary conditions come from the dimensionless 
form of equations ( 32 ) and (33) 1 


D 1 C ° f > G i 
1 . ( 1 0 1 


* N 


i J^o 


K 2 (C°) 3 <C°V - (KjC“) cf' 


= 0 , 


3 1 f9 c - 1 

2. 3( —) ^r?- 

C 2 L 23 J$>0 


f 3 C ; 




= 0 


:=o 


The other two boundary conditions are: 


3. at ^ =1, C 1 =1 , 


for ^ > 0 , 

4. aty«l, C 2 *=C^/cr, forf>0. 


5.2.2 Discretized Equations 


The solution is approximated by a (L+2)th-order Legendre 
polynomial which satisfies the differential equations and 
boundary conditions exactly at N+2 points. Those points are 
chosen to be zeros of the Legendre polynomial of degree over 
the interval (0,1) and with the boundary points 0 and 1. 
Discretizing the spatial derivatives at those points leads to 
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the following set of first-order ordinary differential equations 



i= 2,.N + l. ( 39 ) 



1— 2,#..#N+1, ( irO ) 

where ^ ^, i=2,... N + l are the N inferior collocation points 
while 1 and ^ +2 are exterior collocaxion points which 
correspond to the boundary points ^ =0 and = 1 , 
respectively, and 0= ^ 1 < ^ ^ N+2 = 1. i=2, . .. N + l. 3 j _ , 

is a (N+2) by (N+2) coefficient matrix which comes from the 
discretization of the second derivative ( d 2 C 1 */^'^ 2 or 
O C 2 /S *5 ) at each collocation point. A detailed explanation 
of matrix 3 is shorn in Appendix C. C, (^ .,*£■) and C~ (t ,,t) 
are the concentrations of species 1 and 2 , respectively, ac 
time 'C and position ^ ^. 
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where ^ is an element of another (M-*-2) by (N+2) coefficient 
matrix which comes from the discretization of the first-order 
derivative ( 2>C^ or "^C 2 ) at each collocation points. 

Since only boundary conditions 1 and 2 (equations (43) and (44)) 
require the evaluation of the first-order derivative at ^ 1 =0, 
only the first row of the matrix A, i.e., A, is used here. 

1 I o 

The initial conditions for both cases are the same: 

I.C. 's 

1. at -f=0, c 1 *(^ i .T) = 1 > i = l,...N+2, 

2. at ^ =0, c 2 *('| i ,^)= C^/C^ , i=l,..,N+2. 


5*2.3 Calculation of Current Densit\ 


The current density is related to the flux of the nitrate 
ion at the electrode surface by equation (34). The discretized 
and dimensionless form of equation (34) is; 
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3.3 The Introduction of a Solir.e Point; to the Tiscretited 
Equations 


The diffusion thickness introduced in the derivation 
of the dimensionless form in the previous section can he 
justified by selecting the total electrolysis time such 
that diffusion effects are negligible at x=£ during chat 
period. At the beginning of the electrolysis, the diffusion 
thickness selected according to the total electrolysis time is 
too large. This leads to a concentration profile which is 
fiat in most of the interval with a sharp drop in a very small 
region near *5- =0. This would require a large number of 
collocation points leading to a very stiff set of ordinary 
differential equations. 

This difficulty can be overcome by the use of a method 

f 2 Q \ 

called spline collocation' ' based on the diffusion boundary 
concept. This method maintains a fixed low number of collo¬ 
cation points in the interval which is very small initially 
and will be increased as time increases to account for the 
thickening of the diffusion layer. In the regions outside 
this interval, the concentration is flat. As a matter of 
fact, the spline point can be located as close to the inter- 

( ? Q ) 

face as to achieve any desired accuracy , 

t the concentration gradient 


In order to be sure tha 




at the spline point is zero, the concentration at the Nth 
collocation point, i.e., the last one before the spline 
point, is tested to assure that it is within a very small 
range of the bulk condition. At the time when such a test 
fails, the spline point is moved further into the solution, 

Introducing a spline point,^ s , G< ^ c k 1, a new variable, 
z='V/ "j: s , is introduced. The discretized equations we obtained 
for the case of solution containing cadmium ions are: 


dC 

1 

df la, 


_ ^ r :;+2 * 

* < v 2 > < r* > 21 


i = 2> « e 


a ( >2 ) 

<X+2 * *-i 

T 3,- ,C 2 (a,.T) 

-2K f 

z • 5 s 

^ * ■ ■ - f J ^ J J 

j = l 

l 


^C 2 (z r X)) 


rV o T) 

* 2 


1=2 , .,.h*l. 


For the case of solution containing no cadmium ion, 



c n 
-/ ( 



Equations (43) through (46), which cone from the four ccur.car; 
conditions, become: 

D c b ?l+2 

la ( ~fs } ( )[2I A i fj c i ( s j»T J] + [k 2 (cJ) 3 c 2 ( = 1 ,^) 

- (K-)c; ( t 1>T )] =0, for allt>0. ' (32) 


2. 3( • 


D i . < ?il2 


L >[Z a i,3 C i #(z 3^0 + [E a i,; c 2*^*^ 

2 U j = l ^ c - = l J 


for all X>°* 


3' c i* ( ^ s .t)“ c i* (: :;. 2 't) = for ail T>o, 

4 * 5 2*^ a -T>* s 2 # ‘*a»a»t)*«I/«i. aii t?c 


The initial conditions become: 
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A 


4 


T r ' <3 • 

1 • ^ • O « 

* 

1. at ^=0, C 1 (z^T) 3 1. i=l,...N+2, 

2. at r=0, C 2 *( ZjL , r )= C 2 /C 2 , 1=1, . ..N+2. 


The current density should he changed to s 



5.4 Solution Procedures 


For the case of a solution containing cadmium ions, 
equations (48), (49) and (52) through (55)» with the initial 
conditions, are the working equations used to solve for the 
concentration profiles of N'O^ - and OH - ions. The Cd" 1 "^ ion 
concentration is calculated by the electroneutrality conditic 
in the solution: 
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For the case of a solution containing no cadmium ion, 
equations (50), (51) and (52) through (55) with the initial 
conditions, are the working equations. 

Equation (56) is used to calculate the current density 
at each time / f for both cases. 

5•5 Computer Program Structure 


The computer program for solving the discretized 
equations consists of a main program and 5 main subroutines. 
The first main subroutine evaluates the collocation points 
(roots) of the Jacobi polynomial of order N, as well as the 
firs and second derivatives jf tne polynomial at the roots. 
These derivatives are then used in a second main subroutine 
to caxculate the discretization coefficient matrices A and 3. 
The third main subroutine is used to perform a semi-implicit 


integration using Gear's routine 


(29) 


Inis subroutine cans 


4 external subroutines. The first of these contains the 

explicit expressions (48) and (49) (or (50) and (51) in 

another case) for the discretized coupled first order differentia- 

equations containing the 3. . terms; the second contains the 

-> j 

Jacobian matrix for the same equation; the third performs the 
first stage of Gaussian elimination^ ; of the Jacobian matrix; 
tr.e _ast performs the back substitur.cn of Gaussian elimination 


of 
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of the Jacobian matrix. The first external subroutine called 
by the integration subroutine also calls another subroutine 
coming from the IMSL package. This subroutine solves equations 

(52) and (53) simultaneously to give the surface concentrations 

* * 

of C. and C- . The IMSL subroutine calls another function 

1 c. 

which feeds the equations to be solved. Before writing the 

* * 

function, C ^ (z^,^) (surface concentration of Cp ) is solved 

-g. -fr 

in terms of (z-,^), i = l,.,.Ni-2, and (t-.'T), i = l, . . .11+2 , 

using equation (53) » and then substitute into equation (52) 

* 

to eliminate (z^,^). Thus, only one equation which 

* 

contains only one unknown, is needed. Cnee we 

* * 

obtain (z^,'-£-), C ? (z^,^) can be calculated easily. 

The Jacobian matrices for both cases can also be evalu- 

* 

ated after substituting the expression of (z.,,^) into 
equations (43) and (49) (or (50) and (51) in another case). 

For the case of solution containing no cadmium ion, the Jacobian 
matrix is: 
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which is a 2M by 2N matrix. For the case of solution 
cadmium ions, the Jacobian matrix is the same except 
the terms in the diagonal line of the third quadrant 
subtracted by X and the terms in the diagonal line of 
fourth quadrant are all subracted by 2K(f*(DDXsp/C^ v 


containing 

that 


are all 
the 


T) 


The values obtained from the integration subroutine are 
then used in a fourth main subroutine to calculate the current 
density by equation ( 56 ) . The last main subroutine is used to 
evaluate the concentration at any desired point between 0 and 1 
by interpolation using the values at the collocation points. 

The structure of the whole program is shown in Figure 11. 
Programs from both cases are listed in Appendix D, Programs 
herein were executed on an AMDAHL 4?0/V6 computer. 
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—j i tt; X 
U l *r»i 4.—'-J 

F ITT I?*’G OF THE THEORETICAL 
MODEL WITH EXPERIMENTAL DATA 

6.1 The Determination of the Electrochemical Reaction 
Kinetics Parameters 

The unknown heterogeneous reaction rate cons¬ 
tants wnicn appear in equation ( 52 ) were dexermined 
from xhe currenx-density-vs.-time data obtained in a 
solution which did not contain cadmium ion. The di¬ 
fferential equations ( 50 ) and ( 51 ) which describe the 
transport process were solved by computer to determine 
the rate constants. Comparing equation (31) with equa 
tion (31a), one can see that both K 1 and K 9 depend on 
the applied potential with exponential dependence. 
Thus, X. and ^ should be constant when the applied 
potential is maintained at a constant value during the 
electrolysis. K 1 and K^ were obtained by comparing the 
experimental data in the absence of cadmium ion in the 
electrolyte with the computer-calculated current der.si 

The experimental current-der.sicy-vs ,-time data 


data 



with the electrolyte containing 0 . 005 '^ nitrate ions 
and the applied overpotential of -0,407 was first 
fitted to the computer-simulated data to obtain the 
K 1 and values. Due to the lack of information 
about the values of K 1 and , rough estimates of X. 
and K 0 were evaluated by substituting the applied 
potential V, which is equal to the value of the applied 
overpotential plus the rest potential, to equations 
(31) and (31a). The estimated values of X, and K, were 
increased or decreased until the computer-calculated 
current-density-vs.-time curve had the same shape as 
the same curve obtained experimentally, ('iote, the 
background water decomposition current has been sub¬ 
tracted. ) 

After we obtained the right shape, the anodic 
reaction rate constant, X^, was then fixed. By changing 
X, , one can make the entire curve shift upward. Con¬ 
versely, decreasing the K, value will decrease the 
current density at each time point and thus shift the 
entire curve downward. A set of current-density-vs.- 
time curves with fixed K 7 values and different values 
of X. are plotted in Figure 12. The values of X, and 
K 0 when the concentration of nitrate ions is 0.CC 
and the applied overpotential is -G.^CV was obtained 






Electrolysis Time, ms 
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by interpolating the experimental curve between those 
curves in Figure 12. The values of X 1 and X 2 which 
achieved the best approximation of the experimental 
data are K. = 0.12 X 10~-^ and K^= 0.6 xlO - - 5 , respectively 

These values were then used to predict the cu- 
rrent-density-vs.-time curves at the other ion con¬ 
centrations. Figure 13 shows the computer-simulated 
current density curves and experimental curves at va¬ 
rious nitrate ion concentrations. 

The values of and K 2 at the applied overpo¬ 
tential of -0.60V were different from those at -0.40V. 
The same procedure was used in evaluating the values of 
and K 2 at the condition of -0.60V. Figure 14 shows a 
set of calculated current-density-vs.-time curves with 
fixed K 2 valued and different values of K 1 . The value o 
K„ at the condition of -0,60V was then obtained in the 

i. 

same way as before. The values of K, and K 0 were deter- 
mined to be 0.10 x 10 J and 0 . 50*10 , respectively. 

These values were then used to predict the current-den¬ 
sity-vs . -time curves at the other ion concentrations. 
Figure 15 shews the comparison of the predicted and ex¬ 
perimental current-density-vs,-time curves. 

The reaction rate constants at the applied overpo 
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Figure I'Jt Experimental and Computer-Simulated Curren t-Densi ty-vs . -T ime 

Curve:; with Various Nitrate Ion Bulk Concentrations in the Case 
of FI t-c I rolyte Containing No Cadmium Ion. 
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tential of -0.30V were not pursued. The reason,which 
was discussed previously, is that at this high applied 
overpotential more complicate sequences of reactions 
take place. This simple reaction expression can not 
describe the phenomena. 

The heterogeneous rate expression so determined 
is believed to be a correct one. This claim is supported 
by the fact that the rate constants, and Kg, are 
independent of the bulk concentration at a given poten¬ 
tial. This is indeed the case as was shewn in Figures 
13 and 13 . 

6.2 The Determination of the Homogeneous Precipitation 
Reaction Rate Constant 

Higher current density was observed when cadmium 
nitrate was used instead of potassium nitrate as the 
electrolyte in the electrochemical cell. The hydroxy 
ion produced by the reduction of nitrate ion copreci¬ 
pitated with the cadmium ion in the solution. This 
increased the reduction reaction rate and led to higher 
current. This homogeneous reaction was assumed to be 
linearly dependent or. the degree of supersaturation of 
cadmium hydroxide (see Equation (1^)). Since the beta- 



roger.eous reaction rate constants X„ and were deter¬ 
mined, the homogeneous rate constant d could be evalu¬ 
ated 'ey fitting the second set of experimental data 
with the theoretical model to determine the value of 

Figure 16 shows the current-density-vs,-tine 
curves with various d values while the other parameters 
were held at a constant value. The current-density-vs.- 
time data obtained in a 0 . 0025 M cadmium nitrate solution 
at the applied overpotential of -0.407 is also plotted in 
Figure 16. The reaction rate constant, >, was then eva¬ 
luated approximately from Figure lc. The most probable 
homogeneous reaction rate constant was determined to be 
SO sec -1 . 


k= 100 



ity-vi;.-Time Curves with Fixed 
ua Reaction Rate Coikj tault;. 
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CHAPTER 7 

DISCUSSIONS AND CONCLUSIONS 

7•1 Surface Concentrations of Various Ions as Junctions 
of Time 


The knowledge of the surface ion concentrations can pro 
vide information concerning the electrochemical deposition 
process. Although the surface ion concentrations were not 
available from the experiment, that information could be 
generated by computer simulations. 

Figure 17 shows the surface concentrations of N0^~ and 
0H~ ions as functions of time in the absence of cadmium ion 
in the solution at two potentials, -O.lOV and -0.60V from 
equilibrium potential and at a N0^~ bulk concentration of 
0.005M. In generai,the surface concentration of NO," ions 
decreases as time increases, while the surface concentration 
of 0H~ ions increases as time increases. This is due to the 
electrochemical reaction ( 23 ), in which one nitrate ion 
reacts with two electrons and which produces three hydroxy 
ions. 

For the case of -O.m-CV overpotential, the CH - concen¬ 
tration increases sharply from the bulk value of nearly cere 


KEY 


?4 
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to about 0.5 dimensionless concentration units in less than 1 
msec before the hydroxy ions generated are removed effectively 
by diffusion. Thereafter, the OH - concentration increases at 
a slower rate. For ions, the surface concentration drops 

to about 0.72 dimensionless concentration units from the bulk: 
value of 1.0 in less than 1 msec and then decreases at a 
slower rate when the diffusion can effectively supply the 
reactant for the electrode reaction. The sharp changes of the 
surface concentrations of and CK ions cause the sharp 

decrease of the current density in a short time period in 
the beginning of the electrolysis. 


At a more negative overpotential of -C.60V, the increase 
in surface CH - concentration and the decrease of NO^ concen¬ 
tration are even more significant and last over a longer time 
period. The surface concentration of NC^ drops to about 
0.5 dimensionless concentration units ar.d OH concentration 
increases to about G.c7 dimensionless concentration units in 
2.0 msecs. The behavior is the consequence of the higher 
overall reaction rate at -0.60V. The nitrate ions were con¬ 
sumed at a higher rate which produced more hydroxy ions, 
longer induction time is needed for diffusion to effectively 
supply the reactant from the bulk solution and remove tr.e 
oroduct from the electrode surface. 


Similar profiles of 



suriacs 


cor 




and CH 
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at lower bulk NO^~ concentrations -were also obtained. These 
are shown in Figures 18 and 19* 

When the electrolyte contains cadmium ions, the homo¬ 
geneous precipitation reaction of cadmium ions and hydroxy 
ions near the electrode surface promotes the electrochemical 
reaction on the electrode surface in the cathodic direction. 
This causes the nitrate ions to be consumed at a higher rate 
than the case when no cadmium is present. This is shown in 
Figure 20 where the homogeneous reaction rate constant is 50 

4 

sec~ l . It is interesting to note that the cadmium ion surfac 
concentration increases initially and then starts to decrease 
This phenomenon is explained as follows. Immediately after 
the electrolysis began, N0^~ ions were consumed by the elec¬ 
trochemical reaction which produced an N0^“ concentration 

gradient made the N0^~ ions in the bulk solution diffused 

+2 

toward the electrode surface. The Gd ions, on the other 

hand, moved with the NO^ - ions in the same direction in order 

to maintain the elecroneutrality condition. In the meantime, 

the OH concentration was, however, not so high as to consume 
+2 

all the Cd ions brought in by the transport process. The 
cadmium ions was thus accumulated and resulted in a higher 
value than its bulk condition. As scon as the surface concer, 
tration of CH ion raised to some value, the surface concen¬ 
tration of the Cd ^ icr.s starts ~o decrease because of the 
abundant supply of OH" ions consuming the cadmium ions at 
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this time. 


For the nitrate ions, the surface concentration also 
decreases as time increases at a rate somewhat larger than 
that for the case of solution containing no cadmium ion. This 
is shown in Figure 21. Figure 21 also shows the difference in 
the surface concentration of the hydroxy ion between those two 
cases, i.e., in the presence and in the absence of cadmium 
ions. The initial behaviors of two cases are very similar. 
After this induction period, the ion concentration for the 
case of a solution which contains cadmium ions decreases at 
a more gradual rate as a result of the homogeneous precipi¬ 
tation reaction between cadmium and hydroxy ions. 

7-2 Concentration Profiles of Various Tons at Various Time 


For the case of solution containing no cadmium ion, the 
nitrate ions were consumed and the hydroxy ions were produced 
on the electrode, the concentration gradients were established 
which made the nitrate ions diffused from the bulk solution 
toward the electrode surface and the hydroxy ions diffused 
away from the electrode surface toward the bulk solution. 
Figures 22 through 27 show the concentration profiles of 
and OH ions with different applied overpoter.tials at some 
selected times. From these figures, one can see that the 
diffusion thickness, £ , increases with time. The diffusion 







s 






tn 






CM 






O 




o 


O 




>A 


• 



> 

O 


o 

1 

O 

o 

o 






. 


II 

o 

II 

« 

o 


/-* 


zz 

o 



o 


o 

1 

II 


• 1—1 


•H 


c 


•p 

□ 

-P 

II 

o 




cd 

H 

•H 




u 

cd 

-p 


-P 


-p 

•H 

cd 


C 

* 

c 

«P 

u 


ai 


a; 


■p 


a 


a 

CD 

c 


C 

1 

c 

>P 

0) 

o 

o 

r\ o 

o 

o 

• 

o 

o 

o 

Cm 

C 



z 


i-i 

o 


Pi 


Pi 


a 

II 

rH 


rH 

> 


*T" 

3 


3 

O 


a 

« 

O 

(O 


r-1 





TJ 

3 

rH 

<NJ 



0) 

CO 

cd 

-t- 


+ 

•H 


•rH 

•a 


T3 

r—I 

1 

■P 

o 


U 

A 

<n 

•H 




a 

o 

c 




< 

z 






o r o 



uaui'c '( 7 -HC, 


o o 

_ /_0M„ \ • ,-n r> t" •;» n' ;n' 


4.0 6.0 8.0 10.0 12.0 14.0 16.0 1 

Elec trolysis Time, msec 

21t Comparison of the Surface Concentrations of Nitrate and 
Hydroxy Ions Between Two Oases. 


82 


o 



Os 

o 


co 

o 


N 

o 


sO 

o 


rs 

o 


o 


C\l 

o 


o 

o 


Na 

\ 

x 




Q 





Figure 22 Concentration Profiles of Nitrate and Hydroxy Ions at 
Various Time for the Case of Solution Containing No 
Cadmium Ion. 
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Figure 24s Concentration Profiles of Nitrate and Hydroxy Ions a 
Various Times for the Case of Solution Containing No 
Cadmium Ion. 
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Figure 26» Concentration Profiles of Nitrate and Hydroxy Ions at 
Various Times for the Case of Solution Containing No 
Cadmium Ion. 
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centration Profiles of Nitrate and Hydroxy Ions at 
ious Times for the Case of Solution Containing 




thicknesses at various time are listed in Tables 2 and 3 > 

The concentration profiles of N0^“, OH" and Caions 
at various time for the case of solution containing cadmium 
ions are shown in Figures 28 and 29. The shape of NO ~ 
concentration profiles is the same as that for the case of 
containing no cadmium ion but the N 0 ^ ~ concentration on the 
surface is somewhat lower than that for the case of contain¬ 
ing no cadmium ion. The diffusion thickness increases as time 
increases, 

In the previous discussion of the surface concentration, 
the surface concentration of 0H~ ions increases in the beginn¬ 
ing. During that time, OK" ions diffuse from the electrode 
surface to the bulk solution. Figure 28 shows -chat the 
surface concentration increases and the diffusion layer thick¬ 
ness also increases as time increases. Then the surface 
concentration starts to drop and the diffusion layer thickness 
stays at about the same value. This is shown in Figure 30. 
Table 4 snows the diffusion layer thickness for the diffusion 
of OH" ions as a function of time. In Table 4, one can 
see that the diffusion layer thickness stops increasing as a 
result of the decreasing surface concentration of Oh’" ion. 

+2 

The Cd concentration profiles are shown ir. Figure 29 > 
As the profiles indicate, for each specific time, when the 
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Table 2: Diffusion Thickness as a Function of Time 
for the Case of Solution Containing No 
Cadmium Ion. The Applied Cverpoten'ial is 
-0.40V. The Nitrate-Ion Bulk Concentration 
is Q.005M. 


Time (sec) 


0.1235 X 10 


-5 


0.6121 X10 


-5 


0.1758 X10 


-4 


0.3572 *10 


-4 


0.1033 XlO ~ 3 


0.4737 X 10" 3 


0.1258 X10 


-2 


0.7417 XlO' 


0.3946 XlO 


-1 


0.7032 *10 


-1 


Diffusion Thickness(cm) 


0.25 XlO 


-4 


0.45 X 10" 4 


0.80 X10 


— 


0.15 x10 


0.25 X10 


-3 


0.50 x10" 3 


0.80 XlO' 


0.20 X10 


1 n -2 


0.45 X 10 


-2 


0.60 X10 


-2 
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Table 3: Diffusion Thickness as a Function of Time 
for the Case of Solution Containing No 
Cadmium Ion. The Applied Overpotential Is 
-0.60V. The Nitrate-Ion Bulk Concentration 
Is 0.005M. 


Time (sec) 


0.1313 xio "- 5 


0.4982 X10 


-J5 


0.1642 XIO 


rzr 


0.3318 xio" 


0.1298 xio" 3 


0.4407 XIO" 3 


0.1258 XIO 


-2 


0.8718 XIO 


-2 


0.3229 XIO 


-1 


0.7004 XIO" 


Diffusion Thickness (cm) 


0.25X10 


-4 


0,40 X10 




0.70 X lo -i; 


0.10 X10" 


0.25 xio" 


0.45 XIO " 3 


0.75 xio " 3 


0.25 Xio 


-2 


0.45 XIO' 


0.65 Xio 


-2 
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Cadmium Ion at Various Times 
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30i Concentration Profiles of Hydroxy Ion When the Surface 
Concentration of Hydroxy Ion Starts to Drop. 
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Table 4; Diffusion Thickness of 0K~ Ion as a Function 
of Time for the Case of Solution Containing 
Cadmium Ions. The Applied Overpotential Is 
-0.40V. The Nitrate-Ion Bulk Concentration Is 
0.005M. The Cadmium-Ion Bulk Concentration Is 
0.0025M. 


Time (sec) 

Diffusion Thickness(cm) 

0.1236 XIO' 3 

0.25 X10~ 4 

0.6128 X 10 " 5 

0.45 X10 ~ 4 

0.1798 XIO -4 

0.80 XIO -4 

0.3513 xio " 4 

0.15 xio ' 3 

0.1017 Xio ~ 3 

0.25 XIO -3 

0.4856 X 10' 3 

0.50 X10 ' 3 

0.1148 XIO -2 

0.70 XIO -3 

1 0.8811 XIO -2 

0.85 X 10 ~ J 

0 . 4050 X 10" 1 

0.85 xio ' 3 

0.7031 XlO -1 

0.35 X10 ' 3 
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concentration near the electrode is greater than the bulk 

+2 

concentration of Cd , its value versus distance from the 
electrode decreases slowly in the beginning, then faster, and 
more slowly again before finally reaching its bulk value. At 
the time when the concentration values near the electrode 
surface are still larger than the bulk concentration, the 
concentration decreases to some minimum value at some distance 
from the electrode surface and then increases until it reaches 
the bulk value. The location of the minimum concentration 
point moves toward the electrode surface as time increases. 

At sufficiently long time, the minimum point finally reaches 
the electrode surface. 

7.3 Conclusions 


The transport model for the electrode kinetics and the 
homogeneous reaction presented previously has already explained 
very successfully the characterestics of the current density as 
a function of time for both the case of solution containing 
cadmium ion and that without cadmium ion. This model was then 
used to simulate the concentration profiles for various ions 
involved in the reactions. The behavior of the various ions on 
the electrode surface and in the solution was obtained. The ii 
fusion processes were also known after examining these concen¬ 
tration profiles carefully. These can be summarised as follows 








1. The current density decays with time after a double-layer 
charging. The time for this charging is so short that we 
can not observe it cy the experiment. The decay for she 
current density is a result of the depletion of the nitrate 
ions on the electrode surface. 

2. The current density for the case of a solution containing 
cadmium ions is greater than that for the case of solution 
containing no cadmium ion. This is because the precipitaticr 
of Cd(CH )2 promotes the electrochemical reaction or. the elect 
rode surface by consuming OH” generated in this reaction. 

3. The current density depends on the bulk concentration of 






nitrate ions and the applied overpotential. Increasing 
the NO^ - bulk concentration or the applied overpotential 
will increase the current density. 

k. The applied overpotential is one of the factors which 

changes the reaction rate of the electrochemical reaction 
on the electrode surface. Increasing the applied over¬ 
potential will increase the overall reaction rate. 

5. For the case of solution containing no cadmium ion, the 
surface concentration of N0^" decreases as time increases 
while the surface concentration of OH" increases as time 
increases. This is because NC^“ ions are consumed to produce 
the OH - ions during the electrochemical reaction on the 
electrode surface. 

6. For the case of solution containing cadmium ion, the 

surface concentration of icn also decreases. The 

decreasing rate is, however, faster than that for the c--- 




9 ? 


of solution containing no cadmium ion. The surface 
concentration of OH ion increases in the beginning arc 
reaches a maximum. When the production rate of OH - due to 
the electrochemical reaction is less than the consumption 
rate of 0H“ due to the precipitation reaction, the surface 
concentration of CH _ starts to decrease. To maintain elec- 
troneutrality, the surface concentration of Cd r ^ ions will 
exceed its bulk value and then decrease to the value lower 
than its bulk value. 

7. The diffusion processes for the case of solution containing 
no cadmium ion is that the nitrate ions diffuse from the 
bulk solution toward the electrode surface while the 
hydroxy ions diffuse in the opposite direction. 

8. For the case of solution containing cadmium ions, the 
diffusion directions for both the nitrate ions and hydroxy 
ions are the same as that for the case of solution containing 
no cadmium ion. However, the cadmium ions diffuse toward 
the bulk solution in the beginning and then finally change 
direction toward the electrode surface. 

7 »4 Recommendation far Future Works 

Although the model has successfully explained the cases 
of -0.407 and -G.cQV, iz is not likely valid for the case cf 
-G.oCY due to the unexpected characterestics of the ourrer.t- 
de.nsity-vs.-time curves and the comparison of the current 





density with the limiting case. Vftien the applied cverpoter.tial 
is as high as -0.80V, the reaction mechanisms and electrode 
behavior should be identified by doing some further experiments 
Furthermore, the current-density-vs.-time data at -C.60V applie 
overpotential for the case of solution containing cadmium ions 
is not available yet. It can be obtained by repeating the 
experimental work we decribed before. Once the data is obtained 
the model can be used for comparison with the experimental 
data. 


Finally, the configuration of the electrode in the 
impregnation process is not as simple as we have used in derivi 
the model. It is, in fact, a porous electrode in the real 
case. If a micro-porous electrode is possible to obtain, 
the experiment can be repeated to obtain the current-density- 
vs.-time data. Cnee the data are available, the transport 
processes can be examined for the condition of porous electrode 
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APPENDIX A. 1 


THE "BASIC" MAIN PROGRAM USED IN EXPERIMENTAL WORKS 
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13 

PEN* 



13 

REM* 



20 

REN* 

THIS IS THE BA 

SIC PROGRAM 

— J 

REM* 

TIME TRANSIENT E 

MF'ER: IMENT. 

30 

REM* 

26-0 IS A SAMPLING SUBPOUTIN 

33 

REM* 

LANGUAGE• WHEN " 

MUG" IS EXE 

40 

PEM* 

SAMPLE THE CURRENT OUTPUT A 

- cr 
"tj 

REM* 

INTERNAL. THE ANALOG SIGNAL 

30 

PEM* 

INTO A DIGITAL I 

ATA BY A 12 

33 

REM* 

IN THE MEMORY. V 

HEN THE NI.-N 

60 

REM* 

PEACHES A PRESET 

UALUETHE 

64 

REM* 

TERMINATED. THEE 

E DATA ARE 

6 6 

PEM* 

DECIMAL NUMBERS 

Af iD FR I NTED 

i-.X 

REM* 

THE CAR IMEL2S 

used in thi 

7"0 

REM* 
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REM* 

LOCK S 

ETTING 

c 4 

REM* 

ID DESAMPLING TIME INT- 

“*^T 

REM* 

L,'=i.nnr-iriEL 

NUMBER Ur IE 
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REM* 

SAMP LEZ- 

FROM 

60 

REM * 
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07 POINTS T 

•"2 

REM* 

THE ARGUMENTS 

IN "MUG" SU 

54 

REM* 



•56 

REM* 

IF S=0 THEN 

EXECUTE TEE 

•Z’Zf 

REM* 

IF S=1 THEN 
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'50 

REM* 
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UALU2E PAEE 
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REM* 
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UALUES TO £ 

H -• 

REM* 
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PROGEfiM 


'■0.- 1 15 THE INPUT nrPiiV 
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*54 PEj*1 ; + ; 

58 PEM : * 

57 REM* 

•58 PEM* 

*H'? i^EM******^*****^’** k. r . ...... k.r v , ri- 

170 COfr 4r IG 1 5CRES 4 ' 

18*3 DIM 'A '• 10!) .• V 1 ’! w00*0!.*.. Z■•!3000 •• .• C'•! 500 ■ .■ T 1 300 ' 
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PEM* 
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PEM***** INPUT 
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CFiLL 7HE MUG 
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REM* 



230 

REM-* 1 
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REM* 
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REM* 
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REM* 
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REM* 
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REM 
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340 
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END 
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APPENDIX A.2 

THE ASSEMBLY SUBROUTINE "MUG" USED IN THE MAIN PROGRAM 
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Contents Instructions Coiiiinent: 



0 j l x) C9 RCT Return to main program 
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LISTING OF EXPERIMENTAL DATA 
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:asie 


3.3 Curran*-Density-vs. -Tine Tata for the Case of 
Solution Containing '< r o Cadmium Ion. .titrate- 
Ion 3ulk Concentration* O.QCfM, Initial pH= 3.0, 
Aoplied Overpotential* -0.30V, Rest Potential 
=*-0.35^V. 



OC.S 
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Table 3. 


''•urrer.t-ij er. s i ty - vs 


Solution Containing ' : 


-me -a:a lor tr.e -a 
o Cadmium Ion. :iitr 
Ion Bulk Concentration 5 0.C025M* Initial pH=3.0, 
Applied 0verootential=-0.6C7, Rest Potential 

=*- 0 . 35 ^* 



CiJ 
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m I ’ Pi P 
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■a-ie 3.9 Current-2er.sity-ys ,-Tiae Bata far tr.e rase of 
Solution Containing N'o Cadniun Ion. j’itrat 
Ion Bulk Concentration 3 C.CC125M, Initial 
Applied Overuote.ntiai 3 -O.cCV, Best Potential. 
= -0.354V. 


CUPPEHT I'E’r 15 1T' 
• nf-’P CMC 


. 0001 
. 0002 
. 0003 
. 0004 
. 0003 
. 0006 


. OOOi 7 


. OiOOiO 
. 0009 





104.043 
141.976 
1—3. r?6 
111.177 
102.333 



. 004 

33.3112 

. 006 

3oi. 4433 

. 003 

23.3164 

.01 

22.3116 

.012 

41.0336 

.014 

13.1336 

.016 

13.4043 

. 013 

17.2773 

. 02 

16.3263 

• 0 7, •=; T- 


I 

i 


•*ii fli 
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Table 3.12 Background Current Density of Water. Initial 
pK= 3.0, Applied Overpotential 3 -0,90V, Rest 
Potential 3 -0.35^. 



. 0007 

47.7013 

. OOOS 

43.9433 

. 0009 

41.316 

. 001 

33.6363 

.0011 

37.1344 

.0012 


. 0013 

34.1796 

. 0014 

33.4234 

.0013 

32.3016 

. OO13 

31.3504 

. 0017 

30.7992 

. 0013 

30.043 

. OO 19 

29.6724 

. 002 

23.9212 

.0v34 

22.536 

. 006 

19.9063 

. 003 

13.0233 

.01* 

16.5263 

.012 

15.3996 

.014 

14.6433 

.016 

13.3972 

.013 

13.5216 

. 02 

12.7703 

. 022 

12.3947 

. 024 

11. c.436 

. 026 

11.6436 

. 023 

10.3923 

. 03 

10.3923 

. Uw2 

10.5167 

. 034 

10.5167 

. 036 

9.76561 

. 033 

9.76561 

.04 

9. 39 
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Table 5.1c Currert-Dersity-v; 


- _ me 


z o: 


L3e or 


Solution Containing Cadmium Ions. Nirrare- 
lon Bulk Concencration= 0.C023Mi_Cadni 
Bulk Concentration 5 0.00125M* Initial 

Applied Overpotential = -Q.4CV, Has* Po- 

= -0.35^V. 



■JS P, P 

















O 
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Table 3.19 Current-Density-vs.-T 

Solution Containing Cadmium 


ata tor mhe Case ci 
ns. Nitrate- 
Ion Bulk Concentration* 0. 00125^1 Cadmium-Icr. 
Bulk Concentration* 0.C00o25’fl» Initial pH= 3*Q» 
Apolied Overcotential* -O.^OV, Rest potential 
=‘-0.3 5^ V. 











t-i 













Bulk C on c er. tr a t i on= 0.006251*1, Initial prl= 
Aopiied Overooter.tial= -0.807, Rest Potential 
=*-0.354V. 



M 
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The method of weighted residuals, KWR, is a general 
method for obtaining solutions to equations of change, in 
our case, Fick's second law. In the MWR, one assumes a trial 
function, usually a set of weighted polynomials,substitutes 
this trial function into the differential equation and then 
selects the coefficients of the polynomial terms by specify¬ 
ing that the residual be zero, on the average, at certain 
points. If one evaluates the differential equation at the 
zeros of an orthogonal polynomial, the residual will of 
necessity be exactly zero at these collocation points. 3y 
increasing the number of collocation points, the trial 
function would satisfy the differential equation at more and 
more points^ 1 ^. 



(1) 3.A. Finlayson £ Richard Bellman (Id,), 'mathematics 

in Science and Engineering', Vol. 8?, Academic Press, 
New York, 1972, pp 9. 
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with the boundary conditions of the form 

y(o,t)=g(t) (3) 

y(1,t )=h(t) (4) 

and initial condition such as 

y(x,0)=p(x) (5) 

For a mass transfer problem D would be a dimension¬ 
less diffusion coefficient and f(y) would be homogeneous 
reaction term. Geometric considerations indicate that 

o£ = £3=0 and mathematical considerations indicate that a 

( 2 ) 

suitable but not unique trial function is v ' 

N 

y(x,t)»(l-x)y(0,t)+xy(l,t)+x(l-x)ZIa i (t)P. .(x) (6) 

1 

The ? 4 1 of equation (6) refers to one member of a 
0 ** x 

complete set of polynomials. Table C.l lists some Legendr 
polynomials, the special case of equation (1) -where ct = 3 
and their roots . Note that a Legendre polynomial of degr 
has N real, distinct roots or. the interval 0<x< 1. 


z.k. Finlayson - Richard Bellman (Ed 
in Science and Engineering', Vol, 37, 
New York, 1971, pp 105■ 


'Mathematics 
idenic Press, 
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TA3LE C. 1 


Legendre Polynomial and Their 


KOOtS 


(3) 


fl 

P N -1 

Roots of P.. 

iN ■ 

1 

-l+ 2 x 

0.500000000 

2 

1 - 6 x+ 6 x 2 

0.211324865 

0.788675135 

3 

-l+ 12 x- 30 x 2 + 20 x 3 

0.112701665 

0.500000000 

0.887298335 

4 

1-20x+90x 2 -140x^+70x^ 

0.069431844 

0.330009478 

0.669990522 

0.930568156 


(3) John Villadsen, 'Selected Approximation Method for 
Chemical Engineering Problems', Reproset, Copenhagen, 
Denmark, 1970, pp A3 and A4. 
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Ix is evident that regardless of the value of N cho¬ 
sen one can always rewrite eqn.(6) in the general form 

M + 2 

y(x,t)»y x J " 1 b,(t) (?) 

*- ' ~ O 

<j = l 

One could at this point, as is normally the case in 

using MWR, obtain the coefficients b-(t) by substitution of 

i 

eqn.(?) in the differential equation whose solution is de¬ 
sired. However, the computer programming is greatly simpli¬ 
fied when they are written in terms of the solution of the 
differential equation at the collocation points, y(x.,t) 
where i=i,2,3••.N+2, and each x, is one of the N roots of 
the polynomial or one of the two boundary values. So in 
terms of the N+2 collocation points we may write: 

y(x,,t)= T’xP' 1 c .(t) (8) 

u. , u 

iie now have an approximating polynomial function which 
we will force to fit our partial differential eqn.(2) at 
the N+2 collocation points (i.e. cnoose the b,(t) so that 
our partial differential equation is satisfied). In order 
to do this we will need both the first and second derivatives 
witn respect to x. These derivatives can be calculated ex¬ 
plicitly in terms of x since we have already defined the 
dependence of y upon the spatial coordinate,eqn.(o;and (8). 

It is easily seen that if one considers ail the co- 



13^ 


liocation points, x^, equation (5) represents one center o 
a set of equations: 


N+2 

y(x 1 , i)= y~xj~ 1 'o j( t) (9) 

]=' 



Cne can see that the left-hand side is a vector of 

the trial solution evaluated at the collocation points. 

• « 

Since the term involving xv' - * depends on doth i and a 
matrix will represent this term. The b.(t) will naturally 

O 

be a vector. If vie let a single car represent a vector and 
a double bar a matrix, then the set of earns.(9) through;!! 
can be 'written as one equation in this simplified notation 

y(t/=Q•b(t) lily 

wnere Q.. ; =xi‘~^ 

• u 

The first and second derivatives with respect to x 
may then be obtained from equation ill): 
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— 

dx“ 


(1^) 


where C 


dx“” 1 i , n d 2 x J_1 I 

,d d \*i ° jx. 


it is now si.T.pie to soive lor Ditj tne vector ox o.(t; 

j 

coefficients, since simple rules of matrix algebra apply in 
this notation. Equation (12) yields 




o (t) = Q " y (t) 


Substituting this into eqr.s.(13) and(14) gives 


dy(_ 


dx 


Q' x y(t)= A y(t) 


(15) 


*1 y ( v) _ -\ ,-i ~ 1 ,, ( j. \ _ p rrrzr. 
dx c 


( 16 ) 


where A= C and 3= D Q - ^ or in terms of the i-th collo¬ 
cation point, after some simple matrix multiplication, 


dy(x,t) _ 


MtZ 


dx 


= y : i y V x • , t) 


i J 




, 2 . -. 

— 2^ 1 j J , = T~ 3 --,- c 

.4 > “i,J 


dx 


J =l 
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One can therefore express the spatial derivatives in 
terms of the trial function at the N+2 collocation points . 

Equation (16) can be directly applied tc some diffu¬ 
sion process such as that described by ecns.(2) through (5). 
Equation (2) becomes 


dy (t) _ 
dt 



f(y(t)) 


(19) 


where the vector dy(t)/dt represents the set of time deri¬ 
vatives at the N't 2 collocation points, 'tie note that there 
is no explicit spatial (x) dependence in eqn.(l?). 3y assum¬ 
ing an explicit polynomial in x, eqn.(6), and evaluating it 
at specific values, the collocation points, the problem is 
reduced from a partial differential equation with two inde¬ 
pendent variables to a set of ordinary differential equa¬ 
tions , 


One can expand eqn.(19) into its separate members to 
solve the problem. For the i-th collocation point we can 
write, according to eqr..(19), 



for 


l.C.3. • • 
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One now considers the conditions at the boundaries 
x^=0 and x^^l. Since these values are known from the 
original problem i.e. eqns.(3) and (4) one can subtituxe 
the boundary conditions in eqn.(20) to obtain 


d7 dt’ ~ = *l 3 i.r* (t,+ 3 i,N-2* h(t) ' 




I 3 i,j y(x j» t) ] + f(y(x i’ t)} 

J a| 


for i = 2,3,4.,.N+1. 


One should note that g(t) and h(t) are included so 
that the boundary conditions may be time dependent in which 
case g(t) and h(t) would be known or easily obtained. The 
diffusion problem has been transformed into a sex of N 
simultaneous first order ordinary differential equations 
in N unknowns with initial conditions and can be solved 
numerically by a number of standard methods on a compuxer. 

Evaluation of the current at an electrode, i= nFAD|^-! , 

Sx jx=0 

can be obtained directly from eqn.(17) once the concentra¬ 
tion profile is obtained. 






APPENDIX D 


compute: 


PROGRAMS USING ORTHOGONAL COLLOCATION 


METKCI 






1 

133 



J ********************** ***************«*«**«******4****«*-**««*** 

¥ 'J 


C 




c 

THIS PROGRAM SOLVES EQUATIONS (50) THROUGH (55) n'li.i I »o 

X 


J 

INITIAL COnDI YIGN3 USING SPLINE COLLOCATION M ET h 0 D FG A T :i i 

c 


c 

CASE OF THE ELECTROLYT E CONTAINING NO CADMIUM ION. 

'Z 


- c 

EQUATIONS (DO) AND (51) A HE DISCRZTILED FIRST-ORDER 

C 


c 

ORDINARY DIFFERENTIAL EQUATIONS «HICH ARE SOLVED BY AN 

w 


c 

SEMI-IMPLICIT INTEGRATION SUBROUTINE CALLED DIFSUB USING 

c 


c 

GEAR'S ROUTINE. EQUATIONS (52) THROUGH (55) ARE SOLVED 

'w 



SI J ULI AN ED US L Y 3 Y A SUBROUTINE CALLED l SY3TM Nil 1CH COMES 

c 


V- 

FROM I MS 4> SUBROUTINE PACKAGE u OBTAIN THE SURFACE 

c 


c 

CD NC EN IRA1.0no OF NITRATE AND HYDRO*! ION'S. EQUATIONS (54) 

c 


c 

AND (5 5) ASi In SIR B UL K VALUES. TEE CO.NCENTIi ATION PROFILES 

J 


c 

THUG 03IA1NLD ARE THEN USED TO CALCULATE THE CURRENT 

c 


c 

DENSITY jY : t UA TIG N (5b). TiiE CURRE..T DENSITY DATA ARE 

G. 


c. 

THEN USED TO FIT KITH THE EIP £3Id ENT AL DATA TO IDENTIFY 

J 


J 

THE CATHODIC AND ANODIC REACTION RATE CONSTANTS OF THE 

c 


c* 

SURFACE REACTION. 

<w. 


c 


V. 


c 

- NOdEHCLA.USE - 

- 


N- 

N.(=:.UHBEG oF PoINTS IN X DIRECTION. 



, t 

N=NUM3ER Or' INTERIOR COLLOCATION POINTS. 

c 


r* 

NT = NUMBER OF TOTAL COLLOCATION POINTS INCLUDING TNG 

c 


c 

JCUNDitaY POiNTS. 

n 


v- 

1(1*1) , Y {1 , L) , . . . Y (1 , N) = CON CENT 5 A TIO NS OF NITRATE ION AT 

Z 


r 

S INTERIOR COLLOCATION POINTS. 

z 


v» 

Y (1 , :<♦ 1) ,Y ( 1, & *2) , . . -Y (1,2*N) =CONCEN TR AT 10NS OF H YDa0.vY 

c 


c 

ION AT N INTERIOR COLLOCATION 

c 


u 

POINTS. 

V- 


c 

Cl (1) * Cl (2) cl (NX) CONCENTRATION OF NITRATE ION AT EACH 

c 


c 

POINT IN X DIRECTION. 

c 


s» 

C2 ( 1) , .2 (C) , . ,C 2 (NX) CONCENTRATION OF HYDROXY ION AT tAC.i 

c 


c 

POINT IN I DIRECTION. 

c 


\ z 

X S (1) 3 3U R F.iC o CO NC E N TR AT ICN OF NITRATE UN. 

c 


C 

X3 (2 ) 3 SU R ? AC - CO N"C E N TP AT I ON or HYDaoXY ION. 

J 


! c 

C10= BJLX CONCENTRATION OF NITRATE ION. 

z 


t 

C2 3= 3li LF. C JNC ENTRATI CN OF HYD30CY UN. 

’ w 


c 

D'1AX = 0 IF FUSION THICKNESS, Cd. 

c 


1 

IS 3 LOC AI ION Jr THE SPLINE POINT, 0 < S3 < 1 . 

c 


1 c 

02 3 = i.N CR Ed EN T OF THE SPLINE POINT, EG. 

c 



II MAX 3 MAXUUG .IME TC BE INTEGRATED, SEC. 

c 



JT IdE= Tide. INCaE.’IENT, SEC. 

w 


Z 

T A UMAX 3 D I d E NG Io N LESS MAXIMUM TIME To BE INTEGRATED. 

,'3 


l 

4 rt a* a* O a 1 '1 w 

;; 


c 

» > .\ J — J 1.1 ~ ,i J 2 sj »« L ZZ S I I *1 Z I !IC U Z kj Z r 

■J 


z 

irty iivj.i 3")Z3 Ci’ ."i L D 2 J.( i 10v*• 

J 


1 

G o T A s s r>ACT i o .< u?JE? CF NITRATE ION. 



f ^ 

c o iii j«* i «* o Jo cl u a w«. ro j l\a iu <«j.< oiu a t, x,\.i j jio 



J 

Z *. ii w * 1 G it • 

c 


- 

OO.i 1 = H ol oR D G - -i wO US REACTION RATE CUNGTAnT IN CAT.iJD.C 




Z a *t wC11 w »• • 




CD-CURRENT OE.iSTTY, Ad P/C d* • z . 



-i 

-. • 
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THE DIEFJSIJ.i COEFFICIENTS ■>: ISC NITRATE A.<D HOWH 
IONS ARE 1.<#022-5 A N D 5. 2oOE-5C22/SEJ,.R E3PECTIV EE I. 

i*********************«*’ 


v. 

C 


c 

*c 


IMPLICIT REAL'S (A-H, 0-2) 

DiaENSUN A (20,10) # a {2 □ ,2 o) ,C 1 (oJ) ,C2 (60) ,C (2 0 j , I (6 ,20) , I .'.AX <2 O) , 

1 Dir 1 (20) .DIF2 (2 0) ,DIr 3 (2 0) , ROOT (20) , VECT (20) , ERROR (20) , 

2 S A/2(24,20) , ?N (400) , XS (2) , XX{ 1) , WA l J) ,3A5 (1) 

EXTERNAL A J Xe 

CC.1.1GN/L 1/C J2E 1, C o ZF5 
CO/.;iGN/L2/OOEF2 
CCil.lUN/LJ/N , N r 

co;lio.'«/z.4/a , £ 

COiiaON/L5/\S 
COilEON /Lb/COr.P 4 


C U .El ON/ 7 / .i S. 


no: 


C 0.1A J N / L d / C O E E 3 , C O N 1 , C G N 2 , O A 2 A , 3 S £ A 

cg;«cn/l9/x 


IN P J 4 . DATA -■ 


READ (5,3 10) ND, N, NO, N1,AL, 3E,EPS,TIEA X, IC,DOS 
AEnJ (5,320) C1S,C23,D/iAX,DTj..'lE,Nv,ES 
READ (5 , S J 0) GArfA, B ETA, CCN2, CO.i 1 
B EAD(5,340) j S10,N EE,E P 

-CALCULATE SO.'IE CONSTANTS TO oZ USED IN ONE PSDSRA/.- 

C0EF1 = (1.0 J20-5) /5. 2 b D - 5 
CO.i; 2 = C,. 9/ C 1 £ 

COE?3 = (1. J02D-5) *C1S/D1AX/ES 

COE F4 = 2 . * ( 1.902D-5) *C1 3*9o43 7./0/.AX/ES 

COEr5 = 1./4,ci/ 

D a = D iIA X/ (Ni-1) 

DTAN = ( a. 2 o J - 5) 'DTI XE/0 3 AX/D .1A X 

T A J .1 A X = (5. E o j- 5) *T IX AX / D.l A X/ D.1A X 

SDTAJ=uTA U 

DETA=1. / (NX-1) 

nt=n ♦no*:; i 

NE3=/*N 

-SET INITIAL DI R EN S I ON 4. E SS TI.1Z = 0.-- — 

IAJ= 0. 

-0SZI..E 7 ii L I'AaIITJO A.iJ .11 N l:U 0 D 111 ENS IJ N LES 

. i 4 I N C 4 4fc. N 4 4 0 .4 M J .4 W 0 i N 4 .1 .1 4 .1 4 4 J tl .t 4 .U .1 

SUE joU I. N w D .r j ii5 — — — — — — 

OTAJ .11 = . 0 0 1 * L T A U 
DTrtU.IA =2.0* OTA J 
« ? IT E (u , j j 0 i 


1 





n n r, r, n n n n n n r, n Ci n n r. n n n 
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-E V A LJ AT E THE COLLOCATION POINTS (BOOTS) OF THE 

JACOol POLYNOMIAL OF OLDER N, AS N ELL AS Tci 
FI aST AND SECOND OSHIVAT IV 35 OF THE POLYNOMIAL 
hT TaE HOOTS. - 


CALL JC001lND,N,NO, N1, A L,3 S, DIF1,01F2,DIF3,HOOT) 

WHITE (6,350 j (HOOT (J) , J= 1 , NT) 

-CALCULATE THE DISCA ETISAIION COEFFICIENT MATH 11 A- 

10 = 1 

DO 20 1=1,NT 

CALL DFO? d IND, N , NO ,N 1 ,1, ID, OIF 1, DIF 2, 3173,3001, '/EC I; 

DO 10 J= 1,NT 
10 A (I, J) = /ECT (J) 

20 CONTINUE 

-CALCULATE THE DISCRETISATION COEFFICIENT idlitil 0- 

ID=2 

DO 4 0 I = 1 , NT 

CALL D F 0 ? it (M 0 , S , NO , N 1 ,I,I0,DI?1 , DIF2 , DIF3 , HOOT , V EC. ) 

DO 30 J = 1, NT 
30 3(I,Jj = VSCT (J) 

40 CONTINUE 

-XNpJx THE VALUES OF ARGUMENTS TO 5E USED III THE 

INTEGRATION SUBROUTINE DIF3U3 - 

Mr = 1 

JSTA HI=0 
MAXD S& = 7 
DO 4 5 1=1,N EH 
IMA X (I) =1.0 
45 CONTINUE 
1 = 0 

- INPUT THE INITIAL CONDITIONS- 

DO 5 0 J =1, N 
Y (1 ,J) = 1. 

I (1 , J *N) =CO£F2 
50 CONTINUE 
XS(1) = 1. 

X<j U J *.0.F ** 

WRIISfo, 3*0) 

SO 11 ME = IA U * D .1 A a * D M A X / (5. 2 60-3) 

rtiiD ? R I N . i H r* C J H 3 c* N . D r. N o IT i AT k ;i I j I M r, 
CONCENT 5 A TUNS OF NITRATE AND .i'iZ S 0 X i lo .,5 
COLLOCATION POINTS- 

CALL C J 3 S N T \,C D, A, Y ) 


™ Caii.^.aLA A n 

Pa IN. THE 
A a. .HE NT 
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c 

c 

c 


X 


c 

c 

c 

c 


*RITE(o,3oO) TIHE,CD, X3(1) , X a (2) , K r L A G , IT S R N , a 
Wair£(S,70J) ( Y (1,1) ,! = !,.< EE) 

- iVnxJATS AND PRINT THE COHCZN T RATICNS Or 

nil H AT Z AMD HY D 30 X 'l IONS A 2 EACH POINT 
IN A DIRECTION AT THIS TIHE- 

CALI CONDEN (C1,C2,NX ,N D ,D ET A, ROOT, DIF 1 ,Y) 

WRITE (6,230; (Cl (J) ,0= 1, NX) 

WRITS (6,2-jD) (C2 (J) ,J = 1, NX) 

-C1ECX IS THE INTZSA AilON TI.1S AT THIS .10 A ENT a £ A C n S3 

Itix H A X HI U '1 TINS TO 32 INTEGRATED ——- 

IF (TAU.GE. TAUHAX) GO TO 1 >0 
85 !=!♦ 1 

IF (I. E J. 1) Gu TO 9 0 


x 

C 

C 


IEST IF THE CONCENTRATION G? NITRATE ION AT THE NTH 
COLLOCATION POINT, 12, TH2 LAST ONE SEEDS E THE SPLINE 
POINT IS WITHIN 0.CC01 OxHENSION LESS CONCENTRATION 
JNif OS THE BOUNDARY CONDITION. - 


C 

C 

c 


IF (( 1.-Y ( 1, J) ) . LT. 0. 0001) NO TO 30 

- IF TaS TEST SAILS, THE SPLINE POINT *3 INCREASED, 

HOVED FURTHER INTO THE SOLJTION.- 

CALL CHANGE{CS,DZ5,CO Z S3 , C u Z S C, C J Z S 5 ,FACTOR,D T A J, D T A J HI, 
1 J START) 


D IA J H A , 


EVALUATE THE CONCENTRATIONS OF NITSATs. AND HYDROxY 
10N S AT THc CO L LOC ATI-/ N POI N TS AT T H x N xW COOR DINa T E 
OUx TO THE INCREASE 0? THE SPLINE POINT. - 


C 

CALL EXPAND lY,BOOT,DIF 1,FACTOR,ND) 

C 

c - integrate equations { 5 od and ; 3 i) using the concentsatix 

C Ax Tna wOL xG CATxON ? 0 x N . S AT THIS TIH x TO 0 & T a x »* THE 

C CONCENTRATIONS AT THE SUCCESSIVE II'E - 


9u 


C.ALx D IS S U i. {U ER , TA U, Y, S A VE , DTA J , DT AUU I , 0TA U HA 
x 3 3 0 5 , X F LA G , G S T A 3 T , .1A X 0 E 3 , ? N) 

IS (XSLAG. LT. 1) SC TO 190 
x x rtU 1x=. JO 1 * 0 i A U 
JT A U H A = 2 - 0 * 01A U 
xF(x.Gx.x) jO TO 9 3 

IF U (1,1) .0 1. 1 .. OP . Y (1,2) .or. 1-.3. Y ( 1,3) .01. 


? j , 1 1 r / . H a a , 


.03.ili,x) * , x . i . 







n n n o r i 


1 


1>2 


. OR. OAoS ( T (1,5) - 1. ) . G r. 0- 0304) JU TO 90 
GO TO 93 
9 3 :n=i/ic 
M2=M 1*IC 

IF (I. E .12. OS. TAU. GE.TAUKAX) go TO 95 
GO TO db 

- CALCULATE THE SURFACE CONCENTRATIONS - 

95 XX(1)=Y(1, 1) 

ITSa N= 50 0 

CALL ZSXSTM (A UX2,SP,NS 13, NEC, XX ,11 ESN , • A , B A LER) 

XS (1) 3 XX { 1) 

XS (2) = -i. *CJEF1* (A ( 1, 1) *XX (1) *A( 1, NT) ♦1.0) /A ( 1 , 1} - A ( 1 , NT) /A (1 ,1 ) 
1 *CJcF2 

DO 150 J=1,N 

XS (2) = XS {2) - 3. *CCEF1*A(1,J*1) *Y ( 1, J) / A (1, 1) - A ( 1, J + 1) 

1 «*( 1, J*N)/A{1, 1) 

150 CONTINUE 
GO TO 60 

190 WRIT 2(b, 240) DX,S0TAU, N X, DT AU , DTJ. ME, DM AX , AL, o 5, ZS,EPS,XFLAG, 

1 E?, NSIG 

NnirZ(0,2‘O) SETA, GAMA,CONI ,CJN2 

- FORMATS FOR INPUT AND OUTPUT STATEMENTS- 

2 30 FORMAI(/,26X,*N03-•,1X,9D 1 1.4/30X, 901 1.4/30X, 90 1 1.4/30X, 9 0 1 1.4 

1 /3OX,901 1.4/30X, 90 11.4/3OX, 901 1. 4/30X,*0 11.4/30X, *0 11.4 

2 /JOX,9011.4/3OX,9011.4/3 OX, 9011.4/3 OX, *01 1.4/30 a,90 11.4) 
240 FORMAT (//. 5X, ' X INTERVAL 3 ' ,012.4,32X, ' INITIAL TAU STEP Sj.ZE=' 

1 , 01 2. 4/5X , ' NO OF POINTS IN X DIP* <, 13,2'Jji, 

2 'LaST TAU STEP SIZE 3 ' , 012.4/5 X, ' RE AL TIME INTERVAL 3 ' 

3 ,D10.4,24X, 'DIFUSION LENGTH 3 ' ,010.4//5X,'AL 3 ',F7.2,3OX, 

4 ' u E = ',F7.2,25X,'ZS 3 ',F7.4//5X,'EPS 3 ',D10.2,27X, 

5 ' K F L.kG 3 ' , I3//5X, ' E? 3 ' , 01 0.2,27 X, • NS IG* * , I j) 

245 FORMAT (//,5X, 'REACTION ORDER G? N03- ION =',F7.2,25X, 

1 * a E ACT ION ORDER OF OH- ION = ' , F7.2//5 X, ' 1 ST CONSTANT 3 ' 

2 ,j12.4,25X, '2ND CONSTANT =',012.4) 

2 50 FORMAT(/,25 a,'OH-',21,9011.4/3 OX,9011.4/3CX, 9 0 1 1.4/30X,*0 11.4 

1 /JOa,* Oil. 4/3 OX, 90 11.4/3 OX, 901 1.4/30 X, 9 0 1 1.4/ 33 X , 9 0 1 1. 

2 /3Ua, 9011.4/3 OX,9011.4/3 OX,9011.4/3 OX,*011.4/30X,*011.4) 
310 FORMAT (413,2F5- 1,201 3. 2 ,15, Fd. 4) 

3 20 FORM AT (4012.5,17, F7. 4) 

3 30 FORMAT UF7.r, 2012. 4) 

3 4 3 FORM AT(215,012.4) 

3 50 FORMAT (1 5 X , 1 0 F 1 0 . 6 ) 

3 60 FORMAT (//, 0 1„. 4, 1 X ,0 12.4,2 a. ,2 0 j 0. 1 J 15 , F2 0.7) 

37 0 FoR MAI (/,2 a,oF 2U. 6/2 X, oFJO. oi 
3 90 FORMAT (2X,'CJLLCCATICN POINTS: ') 

3 90 FORM AT l//, 3 a, 'TI9.E ', 12 X ,'CURRENT' ,4bX, ' ION CO NCE NT LATT«.i ' ) 

7 00 FORMAT (3 OX, 1 0 F 1 J . 5/3 Ja, 1JF 1 0. 5) 

STOP 

END 






143 


SUBROUTINE JC03I (N D, N, 30, N i ,AL,32, DIFI , DIE 2, JIF3 , ROOT) 

I M PL IC11 hh X'l* 6 ( A - H , 0- 2 ) 

DIMENSION D iF 1 (ND) ,DIF2 (ND) ,DIF 3 ^ 1«0) , LOOT (ND) 

5U3R0UTINE EVALUATES THE COLLOCATION POINTS (SOOTS) CF THE 
JACO31 POLY NOMIAL OF ORDER X AS NELL AS THE FIRST, SECOND 
AND THIRD DERIVATIVES OF THE POi.Yaj.UAL AT THE ROOTS. 

-IN PUT PARAMETERS- 

INTEGER i'i D= T H E DIMENSION OF VECTORS Di F1 , 3IF2 , DI ? 3 , ROOT. 
INTEGER N = iHE DESS EE OF THE JACOBI POLYNOMIAL, I.E., THE 
NUMBER OF INTERIOR COLLOCATION POINTS. 

INTEGER NO*DECIDES WHETHER X=0 IS INCLUDED AS A COLLOCATION 
POIwT. NO MUST 3 E SET E 2 J A L TO 1 (INCLUDING 1 = 0) 
OR 0 (EXCLUDING THIS POINT). 

INTEGER N1= AS FOR NO, BUT FOR THE POINT 1=1. 

REAL An, 3 2= Til Z PARAMETERS OF THE J AC 031 POLYNOMIAL. 

- OUTPUT PARAMETERS - 

REAL ARRAY EDOT=ONE-DIMENSI ONAL VECTOR CONTAINING ON Ell. 

THE JN-N0 + N1 DEEDS OF THE POLYNOMIAL USED 
IN THE COLLOCATION ROUTINE. 

REAL ARRAY = T HR E E ON E-DIMENSIGNAL VECTORS CONTAINING 

DIF1.3IF2,DIFj ON EXIT THE FIRST, SECOND, AND THIRD 

DERIVATIVES OF THE POLYNOMIAL AT THE EEROS. 


; ******** * * * * **£ 


C -FIRST -VALUATION 0? COEFFICIENTS IN RECURSION FORMULAS - 

C -RECuAjIQN COEFFICIENTS ARE STORED IN DIF1 AND DIF2 - 

C 

A 9= A L* 3 E 
AD=BE-AL 
A?=3 E* AL 

D1F1 (1) = (AD/ (A3*2) * 1) /2 
DIF2 (T)=0 

IF (N ,LT. 2) GO TO 15 
DO 12 1=2,N 
- 1 = 1-1 
-=A3 *2*-1 

DIF 1 (I) = (AS *A D/2/ ( Z* 2) * 1) /2 
IF (I. NE. 2) GO TO 1 0 
DIF2(I) = i»j*rt?*21) /Z/Z /(I *■ 1) 
ju TO 12 
2 = 2*2 

1=21 * (A3* - 1 ) 


10 





c 

c 


lc-ii 


t=i* u * + 'L) 

Oj. 72 (I) s£/„/ii-1) 

12 CO21 INJ E 

--- 30^,1 DOT 27I NA TI 0 N 3? ME* ION MZTHOu Wir.’i 3 'J3 ?F. £3 31 J 2 

OF PLaYIOUSLY DET Z3 7.1 3F J aU0T5 - 

15 X=0 

Du 35 1=1,2 
20 XD=0 
XN=1 
:<D1=0 
X S 1 = 0 

t/O 2 x J = 1 , .i 

X?= (0171 (J; -a) *XN-CIF2 (J) *.c^ 
xpi= (difi (j) -X) *xin-oiF2(j) *xdi-xn 

XD=X2 

xd i=x:m 
X.\ = Xi> 

22 x:n=x?i 

2C= 1 

z*x:;/xa1 

IF (I. E*. 1)30 ro 3 0 
DO 2 5 J = 2,1 

25 uC-ZC-Z/ (X-SOOT (J- 1) ) 

30 L=0/2C 
A = X- - 

IF(DABS(-) .01. 1.D-C9JGC TO 20 
BOOT (I)= X 
X = X».00J 1 

35 continue 

-ADj 2 V 3 21 u A L COLLOCATION POINTS AT X=J 03 1=1 - 

NT = N + NO 4 ;t 1 

.r (NO.iiiy. 1 /} \jk) ro 4 5 
DO 40 1=1,3 

J = .l 4 1 -I 

uo sour (J 4 1 ) =ii0u7 (J) 
nOOT (1) =0 

45 IF (N 1. 2*. 1 ) NOOT (NT ) = 1 

-NO.» 27 A L J AT E DEtiIVATI/Z3 OF POL I 30.11 AL- 

Du 5 0 I=1,3 ? 

(I; 

DIFI (1/ = 1 
0IF2 (1) =J 
JIF3 (X; 

DO o u J = 1 , ,i ! 

* - (j. — • ij j. * 0 3 0 

i~x ~aod r (j) 

DIF3 (Ij = F * 0 11- J (I) 4 3=0172(1, 

DII2 (I, 3 I * J IF 2 (I) «■ 2* DIF 1 (I) 


t 







(i n < i 


14 


3 IF 1 (I) = x *D IF 1(1) 
50 CONTINUE 

fztu a:i 
end 


sJuaouriwZ jfopr (n d, n, no, n i ,i, io, qifi , oifi, dif 3,aooi, v-cr) 
1.1 PL. 10 IT RSAL*d (A-H,0- 3) 

31.IE NS ION 3IF1 (NO) , DIF2 (NO) , DIFS (3 3) ,3uOT (HD) , VSCT (33) 


*******»»»*■ 


C SOjkOUriilS EVALUATES DISCRZTIZA "ION COEFFICIENT MAIJICZS c 

<: AN 3 GA USE j. AN QUADRATURE WEIGHTS, NC.1 ALT ZED TO S J .1 1- C 

c c 

C -INPUT PAR A.'l 2TE S S- C 

c c 

C INTEGER N, JO, N1=AS IN SUBROUTINE CCC31. C 

C INTEGER I = THc INDEX CF THE NODS FOj WHICH THE WEIGHTS ARE v. 

TO ot CALCULATED. C 

INTEGER ID=1N3ICAT0R. 13=1 GIVES THE WEIGHTS FOR DY/3X, C 

13 = 2 GIVES THE WEIGHTS FOR D2Y/JX2, AND 13= J v. 

GIVES THE GAUSSIAN WEIGHTS. THE VALUE OF I IS C 

0 la RELEVANT IN THE LAST CASE. C 

C REAL ASSAY ROOT, C 

C JIF1,3IF2, JlrJ =TH E 0NS - 01 HE NS10NAL VECTORS CONFUTED IN 0 

C SUBROUTINE JC03I. C 

C C 

C - OUTPUT PARAMETERS - C 


C REAL A a 3 AT VECT=THS COMPUTED VECTOR OF «EIGHTS. C 

C 0 

NT=N+N0+N1 

IF (I D. s;. J) Go TO 1 0 
30 20 J = 1 ,J . 

IF (J • NE. I) G0 TO 21 
IF (I 3. nZ. 1) Gu TO 5 
VECT (I) =Di?2 ( I) /DIF1 (I) /2 
Gu TO 20 

5 VSCT (I) =JxFU(I)/DIF1 (I)/3 
GO .o 20 


21 

Y =SO OT (1) 

-LOOT (0) 


VECT (J) =3 

IF 1 (I) /3IF1 (J) /Y 


xr .13- x*. 

2) VZCT(J) = V EC T (J) * (3IF2 (I) /JLF^ 

r > 

CwNliNUE 



jO 4 3 3 0 


10 

l = 0 



O'., 2 5 J = 1 

,i r 


:: = ROui ij) 

A X = X * ( 1 - a ) 

1F{N0.S«. J) A * = \ X / X / X 

(.1 L 4 , . 3) «i X = A \/ ( 1 - X) / (1 - a ) 




1L6 


VEOT (J) =AX/JIF 1 {J) **2 
25 l = Y ♦ / E G T [ J ) 

CO o 0 J = 1, N I 
60 VECI (J) = Vr.CT (J)/Y 
50 HETURN 
END 


sUDRoariNZ intbp (sc,nt,x,soot,di?i,xintp) 

Iii?L 10 IT S2 aL*8 (A-H,0-Z) 

DIMENSION SOoT ( N D) , DI F 1 (NO) ,XINT?(ND) 


-■******, 


««***«*«**», 


c 

C 503SOUTINE EVALUATES THE LAGS AN01 AN I NT 


*******<1***********»****(J 

G 

ESP'JLATION GgEFFIGIINT5 


C - INPUT PASAUSTESS C 

C G 

C INTEGER N I = I tt 2 TOTAL N’KIjES OF COLLOCATION POINTS ( = N+N0 + ..1) C 

G FJa WHICH THE VALUE 0? THE 32P2NDZNT V A SIA RLE Y G 

C IS KNOWN. C 

G HEAL X=TtiZ A3SCISSA X WHERE Y (X) IG DESIRED. G 

g PEAL ASS AY C 

G SOOT,DI?1 =COLLOCATICN POINTS AND DERIVATIVES CF NODE C 

C iJu'i N 0 P. I A L, DERIVED IN GU5RO U TINE JC'j3I. C 

G C 

C - OUTPUT PiaAHETEHS C 

•G C 

G SEAL ASSAY XI NT ?=T H E VECTOR OF INTERPOLATION WEIGHTS. C 

G G 


******* i 


:«#*»♦* *G 


?CL= 1 

DO 5 I = 1 , N T 
Y=X-ROOT (1) 

XINTP (I) =0 

IF lY.E^. 0. 00) XINTP (I) = 1 
5 ?GL= POL * Y 

IF (POL.t*. 0.00) GO TO 10 
DO o 1=1,NT 

o XINTP (I) =POL/DIF 1 (I)/ (X-30oI 11) ) 
10 SETUHN 
END 


jJJnuUi'u 

1 

T.-.p-icir 
j The ns ion 
3 

0 x X g 4* S IJ N 


E OikjJo (*i| » , ’) A / ., d , >i>. x N , *1 ■! .C , up j f ,'i r , i 

JSTA S T, U A X0E A , P N) 

9 u liL ^ 1 [ A * il, ^ “ w ) 

L (0,20) , Y.l AX ( 20 ) , SAVZ(2-*,tJ) , ERROR i E J ) 
A (o i ,P'£PTS7(7,2,j) 

<. UO) ,?JL( 20,20) ,X (2 0) 


-lii 


, 


O J O T 


, ?« k -t 0 J j , 


«T Li rt VJ f 


* * ** ******4rft«*«*att* v ««*a*««*;j 


o n 
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\ 


\ 

i 



i 


c 


S U xl R 0 J Cl HZ xNlZGRAT 
ORDER EQUATIONS AS 
ORE STEP 0? LENGTH 
DECREASED .iliulN TH 
AS LaaGE A STEP AS 
STEP ERROR '# a IC !1 IS 

each component of t 


ES A SET 0? .< ORDINARY DIFFERENTIAL FIRST 
DESCai JED i:i E * U A110 N S 150) AND (51) OVEc 
a at each ca^l. a hay he increased or 

Z A .Ml OS a MI N TO a H A a IN ORDER TO ACHIEVE 
POSSIBLE *SIZE NOT COMMITTING A SINGES 
LARGER THAN EPS IN THE L~2 NOSH, NHxRE 
HE ERROR IS DIVIDED BY THE CO.IRONS l.TS of 


Y MAX. 


THE PROGRAM CALLS FOSS SUBROUTINES NAHSD 
DIFFUN(T,Y,DY) 

H ATiN V (P N, N , M , J) 

SOLVE ( N H , t *'J L, 3, X) 

PEDERV(T,Y,?N,M) 

THE FIRST,DIFFUN,EVALUATES THE DERIVATIVES OF CHE DEPENDiNI 
V ARIAS DCS STORED IN Y(1,I) FOR 1=1 TO N, AND STORES T ;iE 
DERIVATIVES IN THE ARRAY CY. THE SECOND IS CALLED ONLY If 
Hr IS SET TO 1 CS 2 FOR STIFF METHODS. IT ?ESFOR IS A rlaJT 
STAGE OF GAUSSIAN ELIMINATION ROUTINE OF TiiE N 3Y N MATRIX 
STORED IN THE ARRAY PM (M,M) . SOLVE PERFORMS THE jACK 
SIJSSTITUilON PROCESS 0? THE GAUSSIAN ELIMINATION RoUTIN-. 

C u ^ o ‘O ji j u 2 u J 0 il u Y xr Hr To 1, A.ID ^ k x »l £> P a n . x a .« 
DERIVATIVES Of THE DIFFERENTIAE EQUATIONS PSOVI.DED DY DIffJN. 

-PARAMETERS —- 


N= THE NUMBER OF FIRST OR D cS OIFFERsN UAL EQUATIONS. 
1=1115 INDEPENDENT VARIA3L2. 


c 

Y— 

A d 

BY 

N ARRAY CONTAIN 

ING THE DEPENDENT V 

ARIA3LES AND 

c 

W 


THEI 

h 1>Z 

ALSO DERIVATIVE 

S. THE DEPENDENT VARIABLES Ahc. 

c 

z 


S * o £1 

ED I 

N Y (1 ,1) , 1=1 T 

C N. 


V- 


» 

w> A 

V E=A 

DLJ 

CX Of AT L £A ST 

12 *N FLOATING POINT 

^GwAirCIi^ i j w j 

c 

c 


b 

Y TH 

E SUBROUTINE. 



c 


H= 

THE 

3TZ ? 

Si Gw TO 3 E ATT 

EM PIED ON THE NEXT 

STEP. 

w 


HM 

IN=I 

HE MINIMUM STEP SIE 

E THAT MILL BE USED 

FOR THE 


0. 


1 

»i r z j 

RATION. 



c 


:ih 

a a— r 

he a 

A aIM U H SIZE TO 

MHICH THE STEP MILL 

3 w 


c 

E?S = Til 

Z irt.T Jit * Z Z r CCNSTA 

NT- 


c 

J 


= THE 

asr 

HOD INDICATOR. 



Vw 



0 


AN ADAMS PREDICTOR CORRECTOR IS 

USED. 

c 

Z 


1 


A MULTI-STEP 

METHOD SUITABLE FOj 

STIFF S Y S T E M S 


Z 




IS USED. 



c 

Z 


Z 


I HE SAME AS C 

ASE 1, EXCEPT THAT 

T HI S S U 3 R o U T x N _ 

z 

■* 




COMPUTES THE 

PARTIAL DERIVATIVES 

J : xi KUaAw.u 

z 

£ 




DIFFERENCING 

OF THE DERIVATIVES. 


c 

z 

V •/ 

A X = A 

.1 A a \ i ur a L CC AT I 

0N S MuICM CONTAINS 

4. . t X> .! .i l ll j M V> fc 

z 

j 


5 

ACH 

Y SEEN SO FA?. 




z 

EZ 

ROR= 

AN A 

RR A i 0? N ELEME 

.<? 3 

X F.w X 40xlx 


c 



u A Z 

wPEt- ERROR IN Z 

ACH COM PON SN I. 



c 

v ^ 

M .1 xl“ 

A CD 

mp - e r ion code. 

^ xi £, A .J ^ Miiw b » ^ o JOwCMjbr > o • 

z 


J 5 

r 4 a r 

=AN 

l<? r; d i ca tg r 

« 


z 

z 

. . 

.< D ER 

= ruz 

xjAaIMJ'I d e? : V 1 

::vz :.uz 3 .UJZJ 3 z 

V _j - J • 

V. 

c 

PM 

- A d 

uO w«\ 

OF rt T LEAST N* *2 rLOATING puINT LOCAT.ONS. 

z 


1 





************************(2 
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**************** * * ** ***** 


DATA PEHiJi /2.J,4.5,7. 33 j, 10.42, 13.7,17.15,1.0, 

1 2.0, U. 3,24.0,37.39, :>3. 33, 70.0 3, o 7. 37, 

2 3. 0,o.0,9. 1 o7, 12.5,15.3d , 1.0,1.0, 

3 1 2. d,/4.0,3 7. d9, 53-53, 70. 08,37.97, 1.0, 

4 1. , 1. ,0.5,0.1667,0.04133,0.0 0d2t>7,1.0, 

5 1.0,1.0,2.0,1.0,.3157,.07407,.0133/ 

DATA A (2)/- 1.0/ 

IH ET = 1 
Hr LAG = 1 

I? (JSTASI.LE.0) DC TO 140 


1 oo 

DO 

1 10 

1=1,3 





DO 1 

10 J= 1 , K 



1 10 


GA VE (J , I) = 

£ 

(0,1) 


HOLD - 

a 3 2 « 




IF 

(a. 

EQ.HOLD) 

GO 

TO 130 

1 20 

£< acj ft 

= 11/3 OLD 




13 

ET 1 

= 1 




GO 

TO 

750 



1 30 

ilQOL D 

= y w 




T Oi.D = I 

k\ A 0 U *1 — i . 0 

Ir (JSTAdT.GI.O) GC TO 250 
GO TO 170 

140 IF (JSTA.-.T. E*.-1) GC TO 1o0 
-«Q =1 
:i3 = n 

;n = n -* i o 

N2 = Ml + 1 

;; 4 = N ** 2 

:<5 = n 1 

:i o = ii 5 * i 

CALL DIFFJM (T,Y,3A7E (N 2,1) ) 
DO 150 I = 1,0 

1 50 i(2,I) = EA7E(N1+I,1) *.i 

:uiea = J 
K = 2 
GO TO 100 

160 IF ( My. E iJ vO Lu) JSTABT = 1 
T = TOLD 
MO = N^Oi.0 



K - 


♦ i 







GO 

TO 

120 






1 70 

IF 

(IF 

• c* 0 • 


■j 0 

TO 

1 00 



1 F 

( Ou 

- g r • 

=>) 

.j G 

TO 

1 30 



GO 

TO 

(22 1 

■) 

/ «• 

22 , 

223 

,224 , 

22 3, 22 o) ,* 

1 30 

L? 

( 

• J A • 

n 

GC 

TO 

1 30 



AjO 

TO 

(2 11 

0 *• 

12 , 

213 

,214, 

215,216,217) 

1 30 

Kr L 

A j 

= “ *1 







r. zz 

•J3:i 







2 1 1 

A(1) = 

- 1. 

J 






GO 

:o 

2 30 






2 12 

a r- 

) = 


5 J 

J00 

00000 
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A iJ|_ 

= - •). 5 u j 0 00 0 0 0 



JvJ 4 

u 3 3 J 

u 

13 

A (1) 

= - 0 . -t loo 6060600066 b7 



A l 3) 

= -0.750000000 



A (‘0 

- -j. 1 o ob 6 D o 6 6 b 66 b c7 



jO T 

0 Z 30 

2 14 

A (I) 

= -0-375000000 



A (3) 

= - 0-9 1 06 b co 6 6666 6 66 7 



A (4) 

= - 0 . 3333323333 33 3 33 3 



A (5) 

= -0. 0416666666 66 6 66 60-3 7 



30 10 230 

2 15 

A ( 1 . 

= -0.34oo 1 1 1 1 11 1 1 1 11 1 



A iJ) 

= - 1. 0 4 1 66 6 o 5 66 66 7 



A (-) 

= - 0 .4jo 11 1 1 1 1111 1 111 



A p) 

= - o. 10416 6 6066666 66 7 



A (o) 

= -O-OOo J333 3 33 J3 333J3 J 



jO r 

0 2 3 ^ 

1 

1 6 

A ( 1 ) 

= - 0 .32*36 1 1 1 11 1 1 1 1 1 1 



A 13) 

= -1.1 16066660060 66 7 



A 14) 

= -J.o2o000G0C 



A (5) 

= -0.177)3333333333333 



A (o) 

= -J.JziOOOOOCO 




= 0. 00 1 3 3 3 3 6 3 S d d 3 0 3 3 5 0 



Go : 

0 230 

“> 

<- 

17 

All) 

= -u. 3 1o53 13 3 12 16931 2 



A l 3) 

= - 1. ZZ50J000C0 



A (4) 

= -0. 7 5135135 13513516 



A (5) 

= -d.255206333333 33333 



A («) 

= -0.0 436 1111111111111 



A (7) 

= - 0. 0 04 36 1 1 111 11 1 1 1 1 1 11 1 



A (o) 

= -0.0 00 10341 26934 12o634 



oO r J 230 

2 

2 1 

A tl) 

= -I.OOdO'OOOOuO 



inO *. 

0 230 

222 

A (1) 

- - 0 . 000606666660606 7 



A 13) 

= -u. 3j J3J333333333.33 



r* \ •** 
J * 

O 23 ) 

■) 

23 

A ( 1) 

= -d. 54545454 54 54 5-04545 j 



A (3) 

= Ail) 



A (4) 

= - 0 . U0090303G9090 00 *1 



3o r 

0 230 

T 

24 

A ( 1) 

= -O. 40 OOOOOOG 



A (J) 

= - 0 . / OdOOOOOO 



A 1 4 ) 

- - 0 . 2 d 0000000 



A P) 

= -u.J200000JC 



JO 

0 Z J j 


z5 

All) 

= -0. -*37 35 0 20 4 3 7 6 5 c 2 



A (3) 

= - ).3Z)1o7jd22 11o73J 



A 14) 

= -j. 3 lu216 a 731021 J4j 



A ( 5) 

= - 0 . J 5-» 7 44 52 55 -• 7 445 2 c 


— 

-3 •) 

-4 . -t 1 6 3 5J 36 4 36 j 5 J 4 



JO 

10 z Ja 

1 

2 o 

A ( 1 ) 

- -j. -jo 1 6 32553Jc1225 



A (J| 

= -u. * zG 0 3 4 5 2 Go 34 9 Z0 0 
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A (4) = -0. 4 1 o 6606 5 66oo outj o 7 

A (5) = -0- 0 3*20476 19047615 

A (fc>) = -0. 0 1 193476 1904 761 3 

A (7) = -0. 0 00 566 39 34 240 Jo2 3 2 
2 30 X = MQ + 1 
IDuUB = K 
;ity? = (4 - ;i?)/2 
EMQ2 = . 5/FLGAT (S3 ♦ 1) 

BNOi = . 5/F LJ A T(N 0 ♦ 2) 

ZN*1 = .5/T LOAT(UQ) 

PEPBTi - ZP3 

ZU2 = (? Z 31 a i (!J Q, M TY ?, 2) * ? E 23 d) * * 2 
Z = (PE3IST (AC ,«TYF, 1) * PE P'Jii) **2 
ZD2A = (2EoI3T (AO, .HYP, J) *2i?SH) *»2 
IF t 2D«M.2Q.O) GO TC 730 
BN Q = ZP 3*a 30 3/D FLCA I ( N ) 

2 40 I'AZVAL = 111 

go ro (250 /O 30) , ihzt 
2 50 r = r ti 

DO 2 bO J = 2, X 
DO 260 J1 = J,X 

J 2 = X — J 1 * J - 1 
DO 2c 0 I = 1,:,' 

2 60 i (J 2,1) = Y (J 2 ,1) * Y(J2+1,I) 

Do 2 70 I = 1,3 
270 2ac(0a(I) = 0.0 

DO 430 L = 1,3 

CALL DIFFUN (T,Y,SA VS(32,1)) 

IF (IW2VAL. LT. 1) GC TO 350 
IF (3F.L2.2) GO TO 310 
CALL 2EDZ32 (?, Y, FA ,113) 
a = a (i) * d 
DO 230 I = 1, M4 
280 ?»J (I) =22 (I) *a 

290 00 300 I = 1,3 

3 00 2 W (I * (a 3 ♦ 1) - 3 3 ) = 1.0 + ?'.i (I * (Hi* 1)-33) 

I * EV AL = - 1 

CA ul, BAT I a V N , 3 3, J1 , POL) 
i.F (J1.GT.0) GO TO 3 50 
oO TO 440 

310 DO 320 I = 1,:i 

320 SAV^(^,I) = Y(1,I) 

DO 3 4 0 J = 1,a 

a = 22 S* 3.1 AX 1 (EPS, DA 35 (SAV£(4,J) ) ) 



Y ( 1 ,J) = Y (1, J) ♦ ? 
J = A ( 1) *;i/n 
v.A Lw D c c 'J ,1 ( . , (, j A / 
DO 3 j 0 I = 1, a 

z (a 6 , i ) ) 


3 30 

2 * (I «■ (J - 1) * M 3) = 

(3 A 7 £ ( :i 5 ♦ 1, 1 ) - D A V Z (A 1 * I , 1 ) 

/ *■> 

3 40 

Y (I, 0 j = BA 7 2 ( 5, J) 




GO To 2 .>o 



3 50 

it iBF.az.O) go ro 

DO 3 o 0 I = 1 ,!( 

370 


360 

BA;B(*,I) = Y (2, I 

) “ in / u (li I ♦ I , 1) * :i 
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30 

ID -4 

1 J 


3 70 

DO 

3 dO 

I = 

1 , s 

3 80 

JA 7a 

{.45 + 1 

, 1) 

= Y (2, 1) - SAVE {0 UI, 1) *H 


call s 

Da 7a 

(•< 

,?JL,5AVE (N 5+1,1) ,2) 


DO 399 

0=1 , 

.4 


3 99 

SAVE (9 

, 0) ** 

(0) 


4 10 

NT 

= N 




DD 4 20 I = 1,'l 

1(1,1) = i (1.I)+A(1)*3A/S(9,I) 

7(2,1) = T (2,1) - SAVZ{9,I) 

23*02(1) = ERROR (I) + SAVE(9,I) 

1? (DAJS (SA VE (9, I) ) . La. (3:iD*Y.1A.C (I) ) ) t. T = N T- 1 
4 20 03.J TIN 02 

IF (OX.LE.O) GO TO 4*0 
430 CONTINOa 
440 T = TOLD 

If ( (H. La. (1.0000 1 )) . A ND. ( (IWaVAL-.TXYP) . LX.-1) ) uo T 
IF ( (Sc. t w . 0) .03. (I« EVAL. :<£. 0) ) 3ACUM = a ACU A* . 2 5D0 
IWaViiL = :i c 
I3ET1 = 2 
30 TO 7‘jO 
^ 60 KFLAG = -3 
4 70 DO d 0 x = 1 , 

DO 4 d 0 J = 1 , K 
4 60 Y(J,xj = 3A VS (J,I) 

F. = HOLD 
JQ ~ NOOaD 
JSTART =32 
RETURN 

4 90 D = 0- J 

DO 5 00 I = !,.i 

5 00 D = J + (a3FOR (I) / Y M A X (I) ) * * 2 

I * a’/ AL = 0 

IF ID. OT. a) 00 TO 540 
I? (rC.LT.J) 00 TO =20 
DO 5 10 J = 3, K 
DO 5 10 1 ~ 1,N 

5 10 Y {J , i) = Y (J , I} + A (J) *a&303 (I) 

5 20 XFLAO = + 1 
KNEW = II 

IF ( xD J U 3. a a. 1) 30 TO 550 
IDOU3 = ID0J3 - 1 
IF ( IDOUo. 3 X. 1) GO TO 700 
00 5 30 I = 1,11 
530 3* 7E { 1 0,1) = Z3SCF. (I) 

3 0 ID / 0 J 

540 XFLA j = i\: U .lu ~ 2 

IF ( a. wi. (JAIN*1.00001) ) Go TO 740 
T = TOLD 

xr (KFuAu* ui« *5) jC TO 720 
5 50 PP2 = iD/a; **2';i2* 1, 2 
Pa 3 = 1. a * x 0 

x c { [ j J * j u • .1 n a D a E ) mD ft . ( Kc ma^j* yb • ~ 1 j j 'jo TO 3 7 J 
D = 0. 0 


0 4 oO 


1 






152 


DO 5 oO 1= 1 , N 

560 D = J *• ( (2350ft (I) - SA 72 [1 J,I) ) / Yii AX (I) ) ** 2 
?SJ = (D/ZJ?) **EN,}3*1. 4 
570 Pal = 1. 

IF (NQ. L£. 1) GO TO 590 
D = 0.0 

DO 5 80 I = 1, a 

5 80 D = D (7(K,I) /Yii AX (I) ) * *2 
?3 1 = (D/£J«N) **EN Cl *1. 3 

5 90 CONI IN J2 

IF (PP2.L£.P3 3) GO TO 6 30 
IF (P3 3.LT. Pa 1) GO TO 660 

6 00 3 = 1. J/.VUjil (PS 1 , 1. E-4) 

N 2 ,t v = N0"1 
6 10 I DOL'D = 10 

I? ( (KFLAO. £o-1)-AND. (3.LI. (1. 1) J ) GO TO 700 
IF ( NE* *. LE. N J) GO TO 630 
DO o20 I = 1 , N 

6 20 t ( N2* D* 1 , Ij = ZB30S (I) * A ( X) /DFLGAT (K) 
o 30 K - N E ’<i ^ 1 

IF (KFuAG. L*. 1) GO TO 670 
3ACJJ = aAO JO *3 
IBET = 3 
GO TO 7 i 0 

640 1? (NE«j. £*. .J^) GO TO 2 50 
N 0 = N 2 » 2 
GO TO 170 

650 IF (PE2.G1. PL 1) GO TO 600 
NEa'G = 

P. = 1. J/iUAXl (?R 2 , 1. E-4) 

GO TO olO 

6 60 R = 1. J/AilA X 1 (?ti 3 , 1. 3- 4 ) 

N E W v = N ^ ♦ 1 
GO TO 610 
670 13 EX = 2 

P = DM IN 1 (5 , t(8 AX/D A3 8' (H)) 
ii = n«y 

LNErt = H 

IF ( Ni*. N £ * 2) GO TO 6 30 

NO - NSNg 
GO XO 1 7 J 
680 3 1 = 1.0 

DO 6 JO J = 2, K 
21 = 31*6 
DO 6 1 0 I = 1 , “! 

6 90 7 (J »1) = Y (J,I) *31 

IJOU3 = K 

700 DO 7 10 I = J,N 

710 7.1 AX (I j = J.1i.( 1 ( YK AX (I) , D AoJ ( i (1 , I) ) j 

0 OTA 3T = N * 

FZT j NN 

720 17 (N 2 .Ew- 1 1 GO TO 7o0 

CAL- DiFFJit (T,Y,SAVE(M2,1) ) 

3 = H/riO-D 
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oo 7.30 : = i, :i 

i t 1,1) = SAV 12 O , I) 

GAVE(2,I) = :iGLD*5AYE (NCI, 1) 

7 30 i [2,1) = J«Vi(2 , I) *3 

•'4 s 1 
K F LA 0 = 1 
00 I 0 17 0 
7 40 KFLAO = -1 
KNZrf = ii 
J START = Nwi 
SETH ON 

750 TUCGd = OlUXipABS (HHIU/HOO) ,3AC'J.1) 
RAC0H = iMIlil (PAC'JH, DA3S(HHAX/HOLJ) ) 
Si* 1-0 
DO 7 60 J = 2 , X 
a 1 = H 1 *a AC3H 
00 7 o 0 1 = 1 , N 

760 7(0,1) * SAVE (J,I) ■*S 1 

ii = iiOLO-»SACO;i 
60 7 70 I = 1 , N’ 

7 70 1 ( 1,1) = SAVE (1,1) 

I DC 0 0 = i\ 

00 TO (100,250,640),10ET1 
7 30 KFLAG - - 4 
GO TO 47o 

*j»t 0 


SJBBOUTIHE DIFF'JS (T, 77,077) 

IMPLICIT -i Z AL* 3 (A- H, j- Z) 

DICE NS UN 0 77 (23) , 77 (3,20) , A (20,20) ,3 (20,23) ,C (20) ,.(2(1) ,» A (3) , 
1 ^ A3 (20) 


0 * * * * * * : 


sgiuotinz Called 37 cifsob evaluates the derivatives u? 

D Z ? Eli 0 23 I VARIABLES STORED 1.1 77(1 ,1) ,1=1,2*!., .. N 3 STORES 
THE DERIVATIVES 111 THE A3HA7 077. I IS THE Ill LEP Ell 0 Ell T 


/ A i\ X A 3 U £t a 



1 a T - a N A I* A J X 1 

CO;iKwN/u1/COuF1 # COEF5 

Cjr..10:i/L2/J0272 

w 3 /«kl 0N/ u 3 / .i t N T 

C2.1.1 ON /L4/A , o 

C J.ifcl JN/ u 7 / i» j*j i !l 2 . - * »l ? 

ZL = 2 *:i 
J 3 0 J. = 1 / 1 I 

30 v. j »n ^ N J Ci 


0 A mw j 0a r 0 


H - /A w Jwo ^ N ^ .i j- 3 J «? .* A L — i a * )) 




15^ 

ii(i)=YTli, 1) 

i r zs n= 5 ou 

‘A mU Ll *J i, ~J X » i ^ J .C 1 / H r* , | i l Jf , |4 *\ , ^ .1 U , A, 'It 3 ) 

CT=- J. *CoZ7 1 * (A (1 , 1) *Lk (1 ) + A I 1 , j r) * 1. ; / A (1 , 1} - A ( 1, .,2) /« i 1 , 1; ‘U 
DO *4 o i - 1 ^ d 

C (I j = 3(1+1, 1) *CT + 5 (1+ 1 ,NT) +COZ72 
45 CONTINUE 

JTORi IHE D E FIV A TIV Z J IN AsRiii D'{'{ — - 

30 1 20 J = 1 , 

DVi (J) = j (0 + 1 , 1) *.{.( (1) + 3 (3+ 1 , NT) *1 . 

3-Y (J + N) =0 iJ) 

33 110 K= 1 , :j 

DIi l J ) = 3 11 (J) +3 (J + 1 , X + 1) * i' i ( 1, a ) 

3 YY (J+i<) = 3 YY (0+ N) -3 - *3 (J + 1 , 1) *C0E7 1 *A (1 ,X+ 1) *YY (1, A) / A ( 1, 1) 

1 + (3(J + 1, X+1 )-3 (0 + 1 , l) *A (1 ,X + 1)/A (1 , 1) ) *YY ( 1,.(+.<) 

110 CJ NTIN JO 

DYY (J) =302? 1 *C 0 E ? 5 * D Y Y (J) 

Jii (J + N) = J u (J+N) *C02 F5 
120 CON i IN >1L 
F IT'J SN 
EN 3 


j'J-JROJfiN Z PS3ZR V (I, i, ? N J) 

IMPLICIT SEAL* 3 (A - H, Q- Z) 

DIMENSION Y(3,20) , P W (4 00) , A (2 0,2 J) ,3 (20,20) 


■- * a : 


* * V * * * * ■# * * = * < 


: ** * * * Jt * 


C EJSRoUTINE LAL-23 3Y 31:003 STORES DIE PARTIAL 0 
C LZ NI7ATI/Z3 jF THE 3IFFZ3ENIAL Z2JA1ION3 PPOVIjSJ In C 
C SUBROUTINE Dlf’JON. THE PARTIAL DERIVATIVES (TOE JACosiAJ) C 
C -ZPE EVALUATED IN THE LA5T SECTION IN CHAPTER 5. 

C************* = *******»**»»****^=***=Ji***i!t«S*****#4i**»***!»*ii**«ji 


COMM ON/L1/COE?1,CO EF5 
CO.WON/LJ/N , NT 
CC'IMoN/Ln/A , 3 


3o 

1 0 

i\ X = 1 

, J 






30 

4.C 

11 = 

1 , 






P .» 

l II 

+ (2* 

r£ :i“ 

i) * 

■l) 

=C CL 

7 1 *Cu.j?5 *o (H + 1 ,NZ + 1) 


L •< 

(II 

* (-* 

K “ 

i) * 

N) 

— * C C 

Z f 5 * j . * 0 (11+ 1, 1) *COZr 1*\ 

1 1 r Sr. + 1, /A v ■ . " / 

? N 

(II 

+ <.*.! 

* ;j ♦ 

u* 

<:< 

-2) * 

i) =0. 


j p .» 

( II 

+ 2 *N 

* ♦- 

(2* 

K 


) - i 3 (11 *• 1 ,.<•<+ 1 ; - - \I I + ! , 

1) * A ( 1, a .i + 1, / A ( 


10 Cu .J 2 LN J Z 
.',ET j N 
ZN3 


n (i 








! 
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susaoun.iE DiCo:i? (nn,a,'JL, J2) 

DI.1ENS ID .'i A ^20 ,20) ,UL(2 0, 2 0) , SCALES (2 0 ) , IPS (20) 


: ** «*«* * : 


******** *************.2 


SU330UTINE CALLED BY MA7INV PEfiFCdJS A FI8ST STAGE CF 
GAUSSIAN EL IM I NATION SOUTINE. 


C 


1 

2 

J 

4 

5 


1 0 
1 1 
12 

1 3 

1 4 


1 5 


CO Mil ON IPS 
J 2= 1 
N=NN 

DO 5 1=1, N 
I PS (I) = I 

SO NN ail = 0.0 

DO 2 J = 1 ,N 

UL (I, J) = A (I, J) 

I? (SO* NS;i-a ES <U L (1, J)) ) 1,2,2 

aOWN&il = A3 S (UL(I,J) ) 

CO NT I ;i U E 

IF (B^rt.ia.1) 3,4,3 

S CALLS (1) = 1.0/P OW N KM 
GO TO j 

CALL SING 11) 

J 2 -~ 1 

SCALES (I) = 0. 0 
CONTINUE 
Nttl = N-1 
DU 17 K = 1 , N :i 1 
BIG = 0.0 
DO 11 I = K, N 
I? = IPS (I) 

SIuE = A3S (QL(I?, X) ) * SCALES (I?) 
IF (SIDE-BIG) 11, 1 1, 10 

j . j — a i.u 

1DXPIV = I 

.Oil * IN 0 ** 

IF (BIG) 13,12, 13 
CALL SING(2) 

J2 = - 1 

GO TO 17 

IF (IDIP IV -K) 1 4, 15, 14 

J = IPS (K) 

IPS (K) = IPS (IOXPIV) 

(I JiCPIV) = J 
K P = IPS (K) 

PIVOT = JL(X?,K) 

X Pi = K> 1 
DO 1o I * X?1,M 
IP = IPS (I) 

E.1 = - JL {IP , K) /?! 70 2 
J I* (I , a ) * -S3 
do lo j = x?i,:i 
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JL(I? ( J) = UL(I?,J) i- i.1*'JL (HP, J) 

1b CONTINUE 
17 *wQ NT xN J a 
X ? = I P S i N) 

:? 19 , 18,19 

16 CALL SI NO 12) 

J 2 =- 1 
19 PET’J SN 
END 


SUoKOUTiaS SOLVE (JIN ,?UL,3,X) 

IMPLICIT oSAL*8 (A-H, w -Z) 

DI.iaNSi.OJJ PUL (20,20) , 3 (2 J) ,\ (2 J) , IPS (20) 


'*«**#*****#*, 


C SJJ3&OUTIJIE CALLED 32 DIF3U3 PEPrGPHS THE SECOND STAGE C 
C (THE EACH SUbSrirJTICtl PROCESS) OF T.iS GAUSSIAN C 
C ELI Hit! AI ION SOUTINE. 


i v m * + * m* i 


‘ ** *c 


COMMON IPS 
N = JIN 
N?1 = N +1 
IP * IPS i 1) 

X (1) * B IIP) 

DO 2 I * 2, !i 
I P = i?s U) 

I ;i 1 = i-i 

s u ;i = o.o 

DO 1 J = 1,131 

1 sua = sj.i ♦ pol (ip,j)»x(j) 

2 X(I) = a (IP)-SUJ1 
I? = IPS (ji) 

X(N) = X(.«) /POL (I? ,N) 

DO 4 I SACK = 2, N 
I = NP1-ISACE 
IP = IPS (I) 

I ?1 = 1 *] 

S JM a 0. 0 
DO J J = I?1,N 

3 3UJ1 = S0.1 + PUL (IP, J) *.< (J) 

4 X (I) = (aU)-SUH)/PUL (IP, I) 

PSI'U oN 


SUoSJUIlNE SINS (I*HY) 


'*«***« i 


*«*«***«i 




Sa 33 DJTI.it CA-LED JY D EC C HP INDICATES I HE ESSOsS IJJ 


1 


u <J 





15? 

PERFORMING IJi GAUSSIAN ELI/. I NATION aOJTINI. 


v. 

C 


i**************** 


‘*********************0 


1' FOR M AI (54 .i J H ATR IX WITH Z E HO ROW IN DECOMPOSE. 

12 F 0 AM AT ( 0 4 jj 3 S I NS UL A H SAIM a Id DECOMPOSE. ZErtO DIVIDE Id _>CLVL. 

: Id IKPRU7. MATRIX ID NEARLY G *. Li G U L A ? , 


13 

FOLilAr 

(54 HO 

NO CONVERSE 


sour = 

3 




so ro 

(1*2, 

J) , 

IWH Y 

1 

WRITE 

(dour 

,11) 



SO TO 

10 



2 

WRITE 

(doOI 

,12) 



vjO t o 

10 



3 

WRIT E 

i^uU r 

*13) 


10 

RETURN 





u li 0 





GU3AOU 

TINE 

HAT I 

dV (PW , N 


DiilE NS 

ION ? 

a (400) , A(20 

* ** ****** 

** ** * 

*********** 

3UDAJJ7I 

N E CALLED 

JY DIF 


‘ * * * C 


c 


GAUSSIAN E L i M I a AT I 0 d Li GUT Id E OF THE H ATP IX STORED Id 
THE AREA! ?W. 

******** *8 i 


DO 1 I = 1, d 
do 1 J = 1, :i 

1 A (I, J) = II* {J- 1 )*:;3) 
CALL DEC O M P (N , A, U L * J 1) 

RITURN 
END 


jUwRJU 1I.»L CU 5 EdT ( CD , A , Y) 
IMPLICIT RaAL*3 (A-K, 0- Z) 

DI/.E d5 10 d A (20,20) ,Y (8, 20) , XS (2) 


C 31J3ROU TINE CALCULATES THE CURRENT DZNSITi H'i ScUAIIod 

C (5o) . 

v. *"* 

V* W 

c C3=caa«z;«; jl*.oI;y, a^p/c ;<.«*<: : 

C A- OliC jZ1Ar lu CO £ IC I IN I *. ATS I.C , A* C 

u * = c <j j %* w n j it *i x * o j s a r loLuula Ti points* ^ 


***********c 


n r. n n n n n 
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C o M M *yN /Z .j / .t f ,t I 
CuMHON/..b/ XS 
common/Lo/co^e-. 

CD=-CJ&;4* ( A 11 , 1) *X5 (1) *A ( 1 , .• I) ) 
Du bo J = 1. :i 

CD= Cii-ZOH? 4*A(l,J + 1)*i[ ' 1 , J) 

50 CONTINUE 
IETU Rli 
END 


S U 3N j U TI N E CO N C 7. N (C1 ,C 2 , N X , N D, D ET A , E0 01,01F1, l j 
IMPLICIT A E AL *3 (A- r., 0- 1 ) 

DIMENSION Cl l 60) , C2 (60) ,ROOT (2.)) , OlE 1 (20) , Y ( 1,^0) , 13 va, , XT. Nl P (2U) 




C 

c 


S'JSf.OJ TINE EVALJATES THE CO NC ENT u A "I 0 N AT EACH DESIRED 
POINT IN X DIRECTION EY AN IN IZNPOLAT10N ROUTINE JOINS 
THE KN Oil .i CONCENTRATIONS AT THE COLLOCATION POINTS. 

~ — — .A i .I) N ~“ 

C 1 , C2= CD NC E.i 1 LA TIONS u? NITRATE AND n'YDSOXY IONS , 

RES P ECT I / ELY , IN D I HE NS ION L E~>o UNIT. 

D E T A = X IN T E o i A , DIMENSIONLESS, 0 <C D E T A < 1. 


L 

C 


c 

i«;• 


COMM ON/L2/COEE2 
C0MM0N/L3/N,NT 
COMMON/Lb/Xu 
DO 3 0 J=1,N X 
ETA =DETA« (J - 1) 

C A «. L IN YEP ( .< D , N T , ET A , ROOT, Die 1 , X a N1) 
Cl (J) - X i L TP ( 1) *XS (1 ) ♦XINTP (NT) -1. 

C2 (J) = XiNT? ( 1) *X5 (2) +XIN TP (NT) *CO ET 2 
DO 7 0 .•;= 1 ,.« 

- 1 (J) =C 1 (C) +XINT? (X* 1) * Y (1 , K) 

7 0 Co (J) =C2 io) «-X I NT? (K+ 1) *1(1, h+.i) 

30 Cu N T IN U E 
R E TU s.N 
END 


1 ET A'j M A, J S i A.i T; 

2.1 PC IC IT R E .i w * S (A - H, O- E) 


*> ' J - .* J u • ^ u 1 >1 K. e .-'i ^ . J . Hy> j P u . 1 . e > . '1 1 , yw i A .« 0 . .1 . y 4 >l.«i 

A .».i .‘*e* . e* -3 * . v..! A I L ii.\ .i y. J J . . 0 T H w i y C:> eAJ . y .* ~ 

JpI^Ny P 0 I . ■ i. 
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*«****«#* ******* 


****** ******* -1**********************»********♦**£ 


GA7E:>=CS 

IF (OH. i Z . 0 . u 0 1 . A N D . 2 S . L T » J . 0 1 J ) 
IF (- J. G E. J. U 1 O. A N 3 . Z 5. LT- 0. 10 0) 
Ir (u 5 « G Z • 0 - 1 J 0) 3 u 5= 0 • 05 

ZS = 4.S*0C-. 

FACro»l = wZ/SAVZ3 

CO OF J = C0 if j/f ACTOR 

CjZF 4 = 00 if 4/FACTO 3 

CO ZF 3 = 00 Zf 5/F A CTO 2/? AC TO', 

otaj = d:a j/ i o j . 
itaj/.i*. oj i * j ix ’j 

DI AO .1A = 2 . Q * 0 2 A J 
J 5TA &T = 0 
fe :u a u 
end 


0iji — 0 . 0 005 
3iS=0.005 


203230 TINE EIPANO(Y,ROOT,JiF1, F ACT03, NO) 

ILLICIT i>iAu*d ( A - H, 0- C) 

31H E .< 5 i 3 .< t (o ,20) , ROOT (20) , 31F 1 (20) , Yi AT ? * 2 0 ,• , 0 Y (zC ) , ii (2) 


'********* i 




S 0 H300 TI .4 Z S7 A-OATES THE CO NC E 0 T S ATI 0 00 AT TiiE CGLiCCAIlON 

POINTS IN li ii .(EA CO OF 01X AT 2 COE TO THE CHANGE OF THE 

5?LINE POINT. C 

C 

■******«*****w«*»*»******* *ii***»»*»4**y***<.*»*»**** tf , ; 


COHO ON /L 2/C Gi f 2 
C0.'LiUN/LJ/N , NT 
C0.1H0N/L3/:n 
00 7 0 I = 1 , N 
ST (II = 7 ' 1 , I) 

7 0 0 Y (I ♦ :i) = £ ( 1 ,1 f N) 

DO 100 1=1,:) 

Ri=aOJr (!♦ 1) * F ACT CR 

IF(aE.OT.I.O) GO TO 50 

C.iL L INCR/ (N D, N T , BE,20<3T, 01F 1 ,YINTP) 

7 (1 ,I/=tIL 4-/(1) *XS(1) <■ 11.. T P (.1T) *1.0 
7 (1 ,I*N) = TINT ?( 1) *X3 (2) ♦ Cl ..TP (NT) *CJ2F2 
Ou ", ) r. = 1 , .< 

i ( 1 , . j = ; v 1 # I) f Y IN T? (1<* 1 ) 7 () 

a 0 Y 11,1 *■ N ) = Y {1,1 *■ N) Y I N T ? v a * 1) * J 7 (X * N ) 

Go .u 100 
9 0 Y (1 , i j = 1. J 

7 t 1 ,1 * ii) = C 0 i F 2 
100 CCNCIN.Ji 
Ml URN 

£. .« U 


n o 
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F J NOTION AUil (XX , K,JA3 ) 

I.iOLICIT 5 SAo*6 (A- H, *-Z) 

DI IE NS 13 N ti (/0,2 0) ,5 (20,23; , X X (1) , 2 A3 (20) 


C 

c 


«***************< 


t- 


C A J X1 IS THE NAJZ OF THE FUNCTION CALLED 31 ZSTST.1 IN C 

C THE 5U OKJU11N a DIFFUN TO FU So IS a THE VALUES OF EQUATIONS C 

C (52) AND (53) . C 


:****«*’ 


CG.M JN/Ll/COt? 1 ,CO EF5 
COMi1UN/Ld/COEF2 
COU.ION/LJ/N , NT 
C0ilil0N/L4/A, d 

C3dHJN/Ld/-uEF3,CO N1,CON2,SANA,BETA 

2 A3 (20) =-3. *CCEF1* (A (1 , 1) *XX (1) *A ( 1, NT) * 1.0) /A ( 1 , 1) -a ( 1 , :.T) /A (1 , 1 ) 
1 *C J 2f 2 

DO 5 00 J = 1 , N 

2Aa (20) = j AN (20) - 3. ♦COE? 1* A (1 , 0 + 1) *2 A 3 (J) / A (1 , 1) - A (1, J+ 1) 

1 * <A3 (J ♦ N ) / A {1 , 1) 

5 00 CONTINUE 

IF UAS (20) .LT.3.) CAP(20) =0- 

AUX1 =A (1, 1) * XX (1) * A ( 1, NT) *1 . ► (CON2* (QAi (20) ** 3 A HA) -CON 1 * 

1 (XX (1) * * 3 ET A) J/COEF3 

DO (jOO J=1,j 

600 AUX1 = AUX1 *«t ( 1,J* 1) *QA5(J) 

3ETU3N 

END 


FUNCTION A JX2(XX,K,3 AS) 

IuPLIwIT 3 w Ai..* d (A — ri, 0"“-*) 

DIMENSION A (20,20) ,5(20,20) ,I(d,20) ,XX (1) , L A d ( 1) 


■** ***«*«* . 




C 

C A 0 X 2 IS THE NA.1Z OF THE FUNCTION CALLED o'i ZSYETJ la 

C THE CAIN FLuGBAH TO FU3NI5H T rfZ VALUES OF Z<*UATIjNS 

C (52) AND (5d) - 


C 

C 


;**«*»**ae*«<*«a***’ 


CL .IX ON /L 1/COS? 1 ,CO EF5 
C 0 S.'l 0 A/ L2/C UE F2 
C 3.1.10N/L J/N , NT 
C 0 OH O N / L N / A , u 

Co .1.1 ON /L 3/ - O E? J , CO N 1 ,CO N2 , O AC A , BETA 
CJ...1 JN/ui/i 

3aA ( 1) =-J. *C„EF 1* (A ( 1 , 1 ) *Xa (1, ♦ A ( 1 ,NT) *1.0)/A ( 1 , 1) - A ( 1 , NT) /A i ' , 1 ) 
1 *Co E ?2 

30 5 30 J = 1 , N 
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3AS ( 1) =EA ji i 1) - 3- *COE? 1*A ( 1 . J * 1 ) * l V 1 , J) /.i (1 , 1 i - A v 1 1) 

1 *1 ( 1 * J*’ N) /A (1,7) 

500 cociruus 

ic (OrtS (1) . i.r. 0. ) 3.'Bp)=0. 

A U \2 = A I 1 , 1) * A a (1) * A ( 1, :J I) * 1 . ♦ (CC:i2 * ( 3 Aii (1) * * J/» JA) -CO I. 1 * 

1 U-( (1) »*3ETA) )/CO Z F3 

DO 6 JO J = 1 , ,i 

600 AlU2 = A Ji2 *A ( *'L (1, J) 

BiT’JSa 

E:i o 





lOL 


JIUMItlXMXMOtMttXMOtM **«***********««**************«.J 

C c 

C THIS ? ii 0 o R A M oU L7 ^3 z.Q U A T 10 N 3 { 4 d) , ( 4 ■) ) , AND (52) i r. R 0 J*j ii w 

C (55) «xTa I'tiJ INITIAL CONDITIONS USING SPLINE COLLOCATION C 

C METHOD FOR THE CASE CF THE 2LECTR0L ITT CONTA IN I NO CADMIUM C 

C IONS. IN THIS CASE, THE TWO SUBROUTINES, DIFFUN AND PEDLGV, C 

C CALLED BY THE INTEGRATION SUBROUTINE, DIFSU3, HAVE TC BE C 

C CHANGED. MEANWHILE, THE MAIN P20G3AM AS WELL AS THE CTa.i C 

C SUBROUTINES STAY THE SATIS. THE CONCENTRATION PROFILE Or C 

C THE CADMIUM ION :1AT BE OBTAINED BY THE ELECTRONSJT3ALITY C 

C CONDITIO.,. THE CURRENT DENSITY DATA ARE CALCULATED IN IHe, C 

C SAME WAY BY EQUATION (56). THE HETEROGENEOUS REACTION RATE C 

C CONSTANTS OBTAINED BY LAST P3CG2AM ARE USED TO IDENTIFY 


C THE HD .10 SEN £0 US REACTION RATE CONSTANT BY FITTING THE C 

C CALCULATED CURRENT DENSITY DATA WITH THE £.\? S31 HEN TAL C 

C DATA- C 

C 0 

c - ADDITIONAL NOMENCLATURE - C 


C S M A LLK = ri D M 0 G E N E 0 US REACTION RATS CONSTANT, 1/SEC. C 

C D5IGK-0IMENS1ON LESS HOMOGENEOUS REACTION RATE CONSTANT, 

C DB i 0 K = 3 M<t LLX * D MA X ** 2/5. 2b £-5 . I 

C D KS ?= DIM L N 3 10 N L ESS SOLUBILITY PRODUCT OF CADMIUM HYDROXIDE , C 
C DKSP=2.0a-23/C13**2. C 




C** : 


.*«»********=< r 


IMPLICIT REAL *8 (A-H, Q-Z) 

DIMENSION A (2 0,20) ,3 (20 ,20) ,C 1 (60) ,C2 (bO) , C (20) , Y (o ,20) , 

1 0 IF 1 (20) ,DIF2 (2 0) ,DIF3 (20) , HOOT (20) , 7 ECT (20) , E n T. - a 

2 Y MAX (20) ,SAV E (24,20) ,PW (R00) , X3 (2) , X.< {1) , ( J) , GAR 

EXTE RN AL A J XL 

Co MM ON / L1 / ~ U E F1,CO ZF5 
CC '..MON /L2/Z Ot F 2 
C 'ML*, . j , N. 

C . aM O.t / L B / A , L 
COM.1UN/L5/XS 


COMMON/Lo/C OS ?B 
CO MM ON/L 7 /.i GIG, NEE , Z? 
COMMON/LB/EOEF3,CON1,CO N2,G AMA,BETA 
COMMON/LR/» 

CUMMON/L10/0 3 13 K,D KS? 


INPUT DATA — 


•i aA 0 (5,310) 
R EAD (5,320) 
R .AD i 5,3 3 0) 
2 EAD (5,3 4 u) 


.. D , 1, 0 , N1 . ,u, 3.| . .j.. a , , 3 .4 

V-1 3, CEB, DM AX, DTIME, N i , EG 
^ AMA , J ETA , C 0N2 ,CO., 1 , GM ALLK 
N 3 x vj , N Z. , E ? 


. AL C’J LAT E Go HE C0 .*G T AN T. 


: j 


? R j 0 5 a M ~ 


COEF1 = (1. j 02D-5) /5. 26 0-5 
COE F2=C2o/C 1 3 










n r. ( i ( i 


r — 
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COSrJ= l 1. 0 0<.D-5) *C1 S/DOA X/'DS 

0.32 F4 =2. * i 1. 0 D2D-5) «C 1 3«^o4J7./J 1 -iAX/c J 

0 S4 j X = 3 44 ft L L * D.N .A \ * w* A.*. / 5. 2 o ■/“ 3 
DKS?=(F.D-2 3j /(CIS** J. ) 

DX=D.1Aa/(NI-1) 

3In j- (5. iSi-5) 4 DTI MZ/D JAX/DuAa 
TAL' 1 A.X = (5. 2oD-5) *T IX Ai/D.1 A 2/D .1A .( 

SDTA U=DTAU 
d ta = i - / < :i :v -1 ) 
nt=n *o +n i 

:i s«= 2*4-i 

- 2ZT INITIAL DIMENSION-ESS ri.iS-J. - 

7\U~ 0. 


J^4 Ikl U 4. 4 Z 


maximum and .u;;i:iu.i dimensionllss 

i 1 .X ^ I.«C 3 S ME NT TO : : J S Z D IN T H Z IN TEG ? .A 110 .« 
SUBROUTINE D IF S U3- 


d:.\j.u = . do 1 *dta j 
DTAJ /.A - 2 . 0 * DT A'J 
* a II Z ( o , J b J ) 


- £ y * L u A TI THE COLLOCATION POI.'I IS (aCOTS) OF I.iS 

JACOBI POLYNOMIAL OF OBOES N, AS WELL AS IAS 
Fla ST AND SECOND DERIVATIVES Of ThZ PoLi'NuMIA- 
AI r ri £ FOOTS.- 

- .i L lm J C j a 1 ( 4 ! D 14 O .4 J f 4 < 1 f rt m g j 411 i) I r D -i. r 2 , Dir i f j •> o 4 ) 

<< .j i r s (&, j 5 0 ) (?. 0 01 { j j ,0=1 ,»»T) 


ALOJLATZ TiiZ DISCS ZTIZ ATI O.'l COEFFICIS.il OAT! 


ID 

20 


10=1 

DO <. J 1=1/.. I 

C A - L 0 : 0 F A ^ N 0 , N , NO , N 1 ,1 , ID , 01 ? 1 , 01F 2 * D if .1, SOOT , V -Cl) 

DO ID J = 1 ,N T 
A(I,J) = 7-CT(J) 
wul.Tl .ib z. 

"™ ”” — — ^ j.Co L(l* i 4 4k 44 JiJ.44>4i4 4. ..lJkl DD^i 4 4 IwiJk.l 4 4 il 4k' 4 . k'k 2 — 4 .—- 


30 

*0 


DO 4=1. .u 

Cfl .*4 J l' k' ^ - * -4 f kl f 4.4/ f .4 I f 4. f 4 J f J 1 1' lfJlk*^/j4 4 3 f 4. 

D^y JO u — I / ki . 

D (I, J, = 7ZCT (J) 

j., i n.i j., 




1 . 1 PJ. .V A L .J a S Jr AF-jjkik^.T^ 4 j .j Z JSc. 

1 N 4 _ j 3 \T 1 0 N S J j 30 ■) I I.N C J... jJj — — — 


1 
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HF=1 

0 S T A aT = 0 
ISTD?=0 
HAXD£3=7 
DO 47 I = 1 , it Eel 
47 Y.'lAX (I) =1. J 

1 = 0 
n 

C - INPUT THE INITIAL CONDITIONS - 

c 

DO 50 0=1,11 

7 11 , 0 )= 1 - 

50 7(1 ,J+N) =00072 

XS (1 ) = 1. 

XS (2)= C0Er 2 
WRITE (6,J90 ) 

60 TIS£=TAU*DRAa*DHAX/(5. 2 6D-5) 

v- ——————— wALk-ULATE AMD ?8IMT THE CURu—NT DENSIT L AT To..'! — 

0 PRINT THE CONCENTRATIONS 0? NITRATE AND HYDROXY IONS 

C AT TiiE NT COLLOCATICM POINTS- 


call cus en x (c d, a , y) 

"RITE (6 , ioU) TINZ,CD,XS(1) , XS (2) , X F .u A G, IX S3 N , AS 
NRI IE (6,70 Oj (Y {1 ,L) , L=1 , N2R) 

C - EVALUATE AND PRINT THE CONCENTRATIONS 0? 

C NITRATE AND HYDROXY IONS A l EACH POINT 

C IN X DIRECTION AT THIS TIME- 

C 

CALL EONJEN (C1,C2,NX,ND,DETA,ROOT,DIE 1 , Y) 

WRITS (6,200) (Cl (J) ,J= 1, NX) 

WRITE (o,250) (C2 (J) ,0=1,NX) 

C 

C -CHECK IF THE INTEGRATION TIHE AT THIS HU.12 NT REACHES 

C THE HAXillUa TINE TO 3E INTEGRATED- 

C 

IF (IAU. GE. TAUMAX) GO TO 190 
35 1 = 1+ 1 

IF (I. EQ. 1) GO TO 9 0 
IF (ISTDP. E«i- 1) GO TO 90 
C 


n 

c 

c 


TESr IF THE CONCENT RATI ON 0? NITRATE ION AT THE NTH 
E 0 u LCCA TI0 N POINT,IE, T.iS LAST ONE 3EFORZ THE SPLINE 
PC * a I IS * IT HI N J- OOC 1 Dl.i ZNS I JN LESS C D NCEN IR .I w i< 
UNI I OF THE 30UNDARY CONOITION. - 


C 


C 

C 


I? (( 1. -7 (1 , 5) ) . LT. 0. 0C0 1) GO IG )0 


IF THE TEST FAILS, THE SPLINE POINT L3 I NCI LAS E J , IE, 
IJVi„ FURTHER INTO THE SOLUTION. - 












V. 

C 
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lAuL w ri A N vj i [ uS f Du Sf L.C x. 3 t ■- JuC A | D / £A'.I03,DTAU,D.AOHl,DxA J F. t\ , 

J5TA3T) 


■~ —— nVALJATE T iiE CO MC £..« 13 AT ION J J? N IT ? A 7 E ANO HYOSOXY 

IONS AT THE COLLOCATION POINTS AT THE N2X CCGtiOI.UrX 
Oil 2 TO THE INC 3 EASE 0? THE SPLINE POINT.- 


CALL EXPAN 0 (Y ,SOOT,DIF1,FACTOR,NO) 


-INI EG 3 AT E EQUATIONS ('48) AND (4j) USING THE C J NC E N 13 AT ION: 

AT In2 COLLOCATION POINTS AT THIS TI.iE TO OuTAaN I .:E 
CONCENTS AT IONS AT THE SUCCESSIVE TINE - 


90 CALL DIFSU3 (N 2 It , TA U, Y,S AYE, OTA J,0TAUHI , OTA UNA , EPS, NT, l A A X , 

1 £.3303, KFLAG,JSTA3r,HAX0 S3 , ? A) 

IF(XFLAG.LT.0) GO TO 190 
DTAUMI=.001*0TA0 
OTA UM A = 2.0 *D TA r J 
IF (I. GE. 2) Go TO 93 

IF (Y (1,1) .0 2. 1..0R.Y (1 ,2) - IT. 1..03.Y ( 1, 3) . GT. 1. . O.i. i (1 ,4) . ST. 1 . 

1 .03. DAoS ( Y ( 1,5) - 1. ) . GT. 0. 0004) GO TO 90 

GO TO 05 
93 .11 = 1/10 
22=ai*IC 

IF (I.EQ..12.03.TAU.GE.TA 0.1 A X) GO TO 95 
GO TO d5 

-CALCULATE THE S'J? FACE CONCENTRATIONS- 

95 OLDXS2 = XS(2) 

XX(1)=Y(1, 1) 

ITL3 N=500 

CALL ESYSTS (aUX2,EP,NSIG,NEZ,XX,ITE3N,-A,3A3,IER) 

XS(1) =XX (1) 

XS(2) =-i. *COE?1* (A (1, 1) *XX ( 1) 4 A ( 1, NT) * 1.0) /A i 1 , 1) - A ( 1 , N T) /A (1 ,1) 
1 *C0EF2 

00 150 J = 1 , N 

XS (2) =XS (A) - 3. *CCEF1 *A( 1, J* 1) *Y ( 1 , J) / A (1 , 1) - A (1 , J ♦ 1) 

1 *i (1, J«-N)/A(1, 1) 

150 CONTINUE 

IF (XS (2) . LI.OLDX52) ISTO?=1 
GO TO oO 

190 X3IT E (a , 2 40 ) D(,3 OTAU, NX, OTAJ, DTI HE, OH AX, AL,SZ, ES, EPS, .< F EG j, 

1 E?,N SIS 

. SET A, GA M A , CO N 1 , C J N 2, S .1A E L \ 

-F 03 .UTS FOR INPUT A NO UJYPJT ST A 2 El EN TS- 

230 FO AN AI (/,2 5 X, * NO 3- ', IX, :.D 11.4/3 OX, 90 1 1 . 4/3 OX, sJ 1 1. 4 /j OX, 30 11 
1 / iO a , 90 1 1 . 4/3 OX , 9D 11 .4/3 OX, 301 1. 4/3 0 X , 9 s i 1. 4 /j.} ,v , 3D 11 . 4 


.i 


* 





166 


2 

2 40 

1 

2 

3 

4 

5 

2 45 

1 

2 

3 

2 50 

1 

2 

3 10~ 
320 
330 

3 40 
3 50 
360 
3 70 
300 
3 00 
7 00 


/JO 2 , 00 11. 4/302, OD 1 1.4/3 02, 30 11.4/302,9 0 11. 4/JO 2, Ju 11.4) 
FOBMAi (//,5X,'X INTE3VAL=' , 012- 4, 32*, ' INITIAL IAU STEP SUE=' 

,0 12.4/52, ' SO Or' POINTS 10 2 013= ',13,2 92, 

•LAST TAT STEP SUE= ' ,01 2. 4/52,'PEAL TiaE INTSB7AL= ' 
,010.4,242, 'DIFJSION LENGTH 3 • ,010. 4// 52, ' A L 3 ',.-7.2,302 

' 0 li = • ,F7. 2,252 ,'ZJ= ' , F7.4//5 X , ' S?3= ',010.2,27 a, 

' K FL AO= ' , I3//52, ' 2? = ' , J 10.2,27 X,' :i 310= ',13) 

POa.i AT (//, 52, 'REACTION 03D2H OF SC3- ION -'.cl. 2,25a, 

' d 2 ACT 10 N QSDE3 OF OH- UN = ' , F 7.2//S a , ' 1ST CONSTANT =' 

, 0 12.4,252, '2ND CONSTANT =',012.4 
//52 ,'32ACTION RATS CONSTANT 3 ',012.4) 

FOaa AT(/,25 2,'OH-' ,2X, 901 1. 4/30 a,901 1.4/302,001 1.4/302,901 1.4 
/JO a, OD 1 1. 4/3 02 , SO 11.4/3 0 2, 901 1.4/3 0 2,9 011. 4/ jO 2 , 9 0 11 . 4 
/3 02 , 901 1. 4/3 OX, 90 1 1 . 4/3 02, 9 01 1.4/33 2. 93 1 1.4/30 2, ih 11.4) 
FOSHAT {413,2/5. 1, 2010. 2,1 5,F3.4) 


FOR:1AT(401 1.4,I7,F7.4) 

FORMAT k zF7. 2,3012. 4) 

FORMAT (215,012.4) 

FORMAT(152,1 OF 10.6) 

FORM AT (//, j 12.4,1 2 ,D 12. 4, 22,2030. 10,215, F20. 7) 
FORM AT(/, 2 2,qF 20. 6/2 X,6:20.o) 

FORM AT (2 2, ' COLLOCATION POINTS: ') 

FOR a AT (//, j A. 'TI.7S ' , 12 A ,'C0.Ha2NT' ,45 X, ' ION CON 
FCSM AT (302, 10 F 10 . 5/30X, 10F10. 5) 

STOP 


ENTSAllON ' ) 


END 


5 U330U TINS 01FFON (T, YY , DY Y) 

IMPLICIT RSAi.*>3 (A-H.Q-E) 

01 a£ NS ION 0 YY (20) , YY (B, 20) , A (20,2 0) ,3(20,20), C (20) , 22(1) , «A (3) , 

1 ^An (20) 


C SUBROUTINE TALL 2D 3Y 0IF303 EVALUATES THE DERIVATIVES OF 
C DEPENDENT VARIABLES ST 03 EO IN YY (1,1) ,1=1,2*11, AMO STORES 
C THE DERIVATIVES IN THE ASSAY DU. T IS THE INDEPENDENT 
C VARIABLE. 



EXTERN «*L A ; J 21 
COaaON/L1/COEF1,CO EF5 
CuaaJN/L2/COEF2 
coaaoN/LJ/j ,n t 
c j aa ON / L4 / a , 3 
COa.I ON /u 7/.« S. j , N r. E , o? 
CoaaoN/LIO/ObIGK,OKS ? 
I1=2*H 

JO 3 0 1= 1 ,11 
* AS (I) = Y Y (1,1) 

30 CONIINJE 







c 
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CALCULATE THE 7ALJZS 


J U .1 4 IlC w 


V A=0 > 


04 


.4l . 


4‘Yl1) = YY(1, 1) 

IT E B 4 = 3 u 0 

v. t\ L u uo / j i ik (A U a 1 , -- ? / a S a. kj , aA,^Ab,^lb) 

CT=~3. *CGEF 1* (A (1 , 1) *XX (1) *A (1 , .i r) *1. )/A 0 ,1) - A l 1 # NT) /A i 1 , 1) *CCEF2 
DO 43 1=1,4 

C(I) = u (I * 1 , 1) *CT ■*• B { !•*• 1 ,41) *C0£F2 
45 CONTINUE 

-SUEZ THE DERIVATIVES IA A A 4 A Y DYY- 

DO 120 J = 1 , 4 

DYY (J) =o ( J * 1,1) * X i (1) *3( J+ 1 , .«Tj * 1 . 

DYY (J*4) =0 (J) 

DO 110 K = 1 , ;; 

D YY (J) = DYY (J) +3 (J*1, X* 1) *YY ( 1, K) 

110 D YY (J ♦ N > =0 YY (J* 4) -3 . *3 {J *• 1 , 1) *CG£F 1 *A ( 1 ,K* 1) * Y Y ( 1, a ) / A l 1, 1) 

1 * (3 (J * 1, K * 1) - 3 {J ♦ 1, 1) * A (1 , K ♦ 1) / A i 1 , 1 j ) * Y Y ( 1 , «v ♦ 4 j 

DYY (J) =G0£F1 *CCEF 5* DY Y (J ) 

DYY (J «■ 4) =0 YY (J > N ) *C OE F 5 

SO? £3=2. * 0 SIG K * (0.5* (YY {1, J) *Y Y (1, J*- 4) ) -DXSP/Y Y (1 , JO) /Y Y .1 ,0 ♦ 1.) ) 
IF (S li ? E A . 3 1. 0.) DYY (J «■ 4) = D Y Y (J»4) -SUPER. 

120 CJSTIUJc 
it l T 0 b 4 


S036OUT14 L 2ZDE3V (T,7, ?4,43) 

IMPLICIT B E AL*8 (A - H, Q-Z ) 

DIME4S1J4 Y(3,20) ,?W (4 00) ,A (2 0,20) ,3 (2 0,20) 






S 0 330J TI4 £ CALLED 3Y DIPSOS STORES THE PARTIAL 
DERIVATIVES DP THE 0 IF FEW EMI A L STATIONS PbOVIDED I., 

S 'J2HOU TI4 £ THE PARTIAL DERIVATIVES (THE JAC0SLA4) 

4ZF.E EVALUATED 14 THE LAST SECT104 14 CSnPTE? 5. 


C 

C * * 1 


U 

c 

c 


k « * * * »* i 


! ** * * * * = k 


COMM OU/L1 /C02F 1 ,CO EF5 
COMM 04/L J/H , NT 

Ckjktu Ukl / .■* / kk , 3 

C CM.1 04 /l 1 0/0 a IIX , D XS ? 

Du Is rvK= 1 , A 
DO 20 II=1,j 

J* (II + (2 * ab- 2) * 4) =CCEF 1*CC£F5*d (II «■ i , aK ♦ 1) 

P » (II ♦ l 2 * E a- 1) * 4 ) =-CC£?5*J . * 3 (11 ♦ 1 , 1) *C 0 E r 1 * A ( 1 , X;. * 1) / A k ' , 1 ; 

IF l II. Z*. K£) GO TC 30 

Jkj k D 4y 

30 Sl?=2.*JjT0K*(G.5*(Y(1,II) ♦ Y ( 1 , 11 ♦ 4 ) ) - D .\«, ?/ Y v 1 , 11 ♦ 4 j / i k 1 , j 

IF(JU?.UT.O.) P * (II ♦ ( 2 *<K- 1) ) = 2 * (11 *■ ( 2 * £ K - 1) * :.) - OS I G v . 

PA ( II*2»a*4* (2 * K a - 2 ) *4) =0. 



40 





1 
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?* (II 0 * H * »< ♦ (2* XK — 1)*M)-(.JiII*1,i\.;*'1)-b(II*1 / 1) *A(1,Ki\+1)/Al1 f 1)) 

♦C0SF5 

Ir(II.i^.KjC) JO TC u 0 
JO TO 20 

60 IF( jjp. ji. o.) ?w (n*2*M*.Sf {2*K\- ij *:i) =?* u*.\K-1) * jj -2. * 

1 0313 X * (. 6 *-2. *0ao ?/T ( 1, II +:•) / 

2 i (i, 11 + :i) / y (i, 11 *• j ) j 

20 CUSTI.iJJ 
10 CO.IT IN US 
l\ 2 T J d N 
END 


***************,, 


»******♦*«(; 


THE CTrlEii J J 3 30 'IT Itf E3 CAL1.20 J'{ THU PSuJTAH STAY T.iE SA.1E C 

as those ::i the pisst peoosah. c 


*******«*************r 


* 








